text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE PartialTypeSignatures #-}
module Main where

import Test.Tasty
import           Test.Tasty.HUnit
import Data.List
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Storable as VS
import Bio.Data.Bed.Types
import Data.Maybe
import Conduit
import qualified Data.Matrix.Static.Dense as D
import qualified Data.Matrix.Static.Generic as D
import qualified Data.Matrix.Static.Sparse as S
import Data.Matrix.Dynamic (Dynamic(..), matrix)
import Data.Matrix.Static.LinearAlgebra
import Statistics.Sample hiding (covariance)
import System.Random.MWC (create)

import Bio.ChromVAR
import Bio.ChromVAR.Utils

main :: IO ()
main = defaultMain $ testGroup "Main"
    [ testCase "sortBed" test
    , testCase "sortBed" bgTest]
    --, testCase "sortBed" peakTest ]

test :: Assertion
test = (map round' $ concat dev', map round' $ concat z') @=?
    (map round' $ concat deviations, map round' $ concat z)
  where
    [(dev', z')] = runIdentity $ runConduit $ yieldMany [(4, cellByPeak)] .|
        computeDeviation 3 3 expectation bg peakBymotif .| sinkList
    expectation = let e = D.replicate 1 @@ (S.fromTriplet (U.fromList cellByPeak) :: SparseMatrix 4 3 Double) :: Matrix 1 3 Double
                      s = VS.sum $ D.flatten e
                  in D.toList $ D.map (/s) e
    round' :: Double -> Double
    round' x = fromIntegral (round (x * 1e4)) / 1e4

cellByPeak :: [(Int, Int, Double)]
cellByPeak = [(0,1,1), (0,2,2), (1,0,1), (1,1,1), (2,0,2), (2,2,1), (3,1,1), (3,2,1)]

peakBymotif :: [(Int, Int, Double)]
peakBymotif = [(0,0,1), (0,1,1), (1,0,1), (1,1,1), (1,2,1), (2,0,1), (2,2,1)]

peaks :: [BED3]
peaks = map (\(chr, i) -> BED3 chr i (i + 500)) dat
  where
    dat = [ ("chr1", 76585873)
          , ("chr2", 42772928)
          , ("chr2", 100183786) ]

bg :: [[Int]]
bg = [ [1,0,1]
     , [1,0,1]
     , [1,2,0] ]

deviations = transpose [ [0.1728395, -0.4444444, 0.1728395, -0.07407407]
             , [-0.2910053, 0.3174603, 0.2116402, -0.19841270]
             , [0.7407407, -0.6349206, -0.7407407,  0.63492063] ]

z = transpose [ [1.1547005, -1.1547005, 1.1547005, -1.1547005]
    , [-0.5773503, 0.5773503, 0.5773503, -0.5773503]
    , [3.2331615, -1.1547005, -4.0414519, 9.2376043] ]

bgTest :: Assertion
bgTest = do
    input <- readData "tests/data/coordinates.tsv"
    expected <- readData "tests/data/transformed.tsv"
    let actual = case matrix (map (\(a,b) -> [a,b]) $ U.toList input) of
            Dynamic mat@(D.Matrix _) -> U.fromList $ map (\[a,b] -> (a,b)) $
                map VS.toList $ D.toRows $ whiten Cholesky mat
    U.map (\(x,y) -> (round' x, round' y)) actual @=?
        U.map (\(x,y) -> (round' x, round' y)) expected
  where
    readData fl = U.fromList . map (f . words) . lines <$> readFile fl
      where
        f [x,y] = (read x, read y)
    round' :: Double -> Double
    round' x = fromIntegral (round (x * 1e4)) / 1e4

peakTest :: Assertion
peakTest = do
    coord <- readData "tests/data/coordinates.tsv"

    expected <- map (map (\x -> coord U.! (read x - 1)) . words) .
        lines <$> readFile "tests/data/background.tsv"
    actual <- fmap (transpose . map (map (coord U.!) . U.toList)) $
        create >>= getBackgroundPeaks 50 coord

    print $ zipWith f expected actual

    --print $ flip map input' $ \x -> meanVarianceUnb $ U.fromList $ map dist $ comb x
    return undefined
  where
    f xs ys = meanVarianceUnb $ U.fromList $ map dist $ map (\[a,b] -> (a,b)) $ sequence [xs, ys]
    readData fl = U.fromList . map (f . words) . lines <$> readFile fl
      where
        f [x,y] = (read x, read y)
    comb (x:xs) = zip (repeat x) xs ++ comb xs
    comb _ = []
    dist ((x1,y1), (x2,y2)) = sqrt $ (x1 - x2)**2 + (y1-y2)**2

max_stars_repo_name: kaizhang/ChromVAR

max_issues_repo_name: kaizhang/ChromVAR

max_forks_repo_name: kaizhang/ChromVAR


text:
import Data.Complex (Complex((:+)))
import qualified Data.Complex as C (realPart)
import qualified Data.Map as M (fromList, Map)

main = print $ mandelbrot (-2, 2) (-2, 2) 20 20 40

quadratic :: Num a => a -> a -> a
quadratic c z = z^2 + c

isBounded f lim z k
    | lim < k = True
    | C.realPart z <= 2 = isBounded f lim result (k+1)
    | otherwise = False
    where result = f z

isMandelbrot c lim = isBounded (quadratic c) lim 0 0

mandelbrot (x1, x2) (y1, y2) n1 n2 lim = M.fromList l
    where points = [ (x, y) | x <- linspace (x1, x2) n1,
                              y <- linspace (y1, y2) n2 ]
          pointMandel point = isMandelbrot (uncurry (:+) point) lim
          l = [ point | point <- points , pointMandel point ]

linspace (x1, x2) n = x1 : linspaceNext difference [x2] (n - 2)
    where difference = abs (x1 - x2) / (n - 1)

linspaceNext d li n
    | n == 0 = li
    | otherwise = linspaceNext d (ele : li) (n - 1)
    where ele = head li - d

ext: hs

max_stars_repo_name: zhiyanfoo/complex-sets

max_issues_repo_name: zhiyanfoo/complex-sets

max_forks_repo_name: zhiyanfoo/complex-sets


text:
module Gradient
    (numericalDiff
    ) where

import Numeric.LinearAlgebra

numericalDiff :: Fractional a => (a -> a) -> a -> a
numericalDiff f x = (f (x+h) - f (x-h)) / (2*h)
    where h = 1e-4


text:
{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeOperators       #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
module Test.Grenade.Layers.Internal.Transpose where

import           Grenade.Layers.Internal.Transpose

import           Numeric.LinearAlgebra                  hiding (konst,
                                                         uniformSample, (===))

import           Hedgehog

prop_transpose4d_behaves_as_reference_swap_three_axis = withTests 1 $ property $ do
  let input           = (12 >< 4) [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47 ]
      perms           = vector [2, 0, 1, 3]
      dims            = [2, 2, 3, 4]
      expected_output = (12 >< 4) [ 0, 1, 2, 3, 12, 13, 14, 15, 24, 25, 26, 27, 36, 37, 38, 39, 4, 5, 6, 7, 16, 17, 18, 19, 28, 29, 30, 31, 40, 41, 42, 43, 8, 9, 10, 11, 20, 21, 22, 23, 32, 33, 34, 35, 44, 45, 46, 47 ]

      output = transpose4d dims perms input

  output === expected_output

prop_same_pad_pool_behaves_correctly_at_edges_three_channels = withTests 1 $ property $ do
  let input           = (12 >< 4) [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47 ]
      perms           = vector [3, 2, 1, 0]
      dims            = [2, 2, 3, 4]
      expected_output = (24 >< 2) [ 0, 24, 12, 36, 4, 28, 16, 40, 8, 32, 20, 44, 1, 25, 13, 37, 5, 29, 17, 41, 9, 33, 21, 45, 2, 26, 14, 38, 6, 30, 18, 42, 10, 34, 22, 46, 3, 27, 15, 39, 7, 31, 19, 43, 11, 35, 23, 47 ]

      output = transpose4d dims perms input

  output === expected_output

tests :: IO Bool
tests = checkParallel $$(discover)

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE FlexibleInstances #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Foreign.VMath.VFloating
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Vector Floating operations.
--

module Foreign.VMath.VFloating (
    VFloating(..)
    ) where

import Foreign( Ptr, Storable, peek, poke, advancePtr )
import Data.Complex( Complex(..) )

import Foreign.VMath.VFractional
import Foreign.VMath.Double
import Foreign.VMath.Zomplex

-- | Types with vectorized 'Floating' operations.
class (VFractional a, Floating a) => VFloating a where
    vExp :: Int -> Ptr a -> Ptr a -> IO ()
    vSqrt :: Int -> Ptr a -> Ptr a -> IO ()
    vLog :: Int -> Ptr a -> Ptr a -> IO ()
    vPow :: Int -> Ptr a -> Ptr a -> Ptr a -> IO ()
    vSin :: Int -> Ptr a -> Ptr a -> IO ()
    vCos :: Int -> Ptr a -> Ptr a -> IO ()
    vTan :: Int -> Ptr a -> Ptr a -> IO ()
    vASin :: Int -> Ptr a -> Ptr a -> IO ()
    vACos :: Int -> Ptr a -> Ptr a -> IO ()
    vATan :: Int -> Ptr a -> Ptr a -> IO ()
    vSinh :: Int -> Ptr a -> Ptr a -> IO ()
    vCosh :: Int -> Ptr a -> Ptr a -> IO ()
    vTanh :: Int -> Ptr a -> Ptr a -> IO ()
    vASinh :: Int -> Ptr a -> Ptr a -> IO ()
    vACosh :: Int -> Ptr a -> Ptr a -> IO ()
    vATanh :: Int -> Ptr a -> Ptr a -> IO ()

    vExp = vop exp
    vSqrt = vop sqrt
    vLog = vop log
    vPow = vop2 (**)
    vSin = vop sin
    vCos = vop cos
    vTan = vop tan
    vASin = vop asin
    vACos = vop acos
    vATan = vop atan
    vSinh = vop sinh
    vCosh = vop cosh
    vTanh = vop tanh
    vASinh = vop asinh
    vACosh = vop acosh
    vATanh = vop atanh

vop :: (Storable a, Storable b)
    => (a -> b) -> Int -> Ptr a -> Ptr b -> IO ()
vop f n src dst | n <= 0    = return ()
                | otherwise = do
    a <- peek src
    poke dst $ f a
    vop f (n-1) (src `advancePtr` 1) (dst `advancePtr` 1)
{-# INLINE vop #-}

vop2 :: (Storable a1, Storable a2, Storable b)
     => (a1 -> a2 -> b) -> Int -> Ptr a1 -> Ptr a2 -> Ptr b -> IO ()
vop2 f n src1 src2 dst | n <= 0    = return ()
                       | otherwise = do
    a1 <- peek src1
    a2 <- peek src2
    poke dst $ f a1 a2
    vop2 f (n-1) (src1 `advancePtr` 1) (src2 `advancePtr` 1)
         (dst `advancePtr` 1)
{-# INLINE vop2 #-}

instance VFloating Double where
    vExp = vdExp
    {-# INLINE vExp #-}
    vSqrt = vdSqrt
    {-# INLINE vSqrt #-}
    vLog = vdLog
    {-# INLINE vLog #-}
    vPow = vdPow
    {-# INLINE vPow #-}
    vSin = vdSin
    {-# INLINE vSin #-}
    vCos = vdCos
    {-# INLINE vCos #-}
    vTan = vdTan
    {-# INLINE vTan #-}
    vASin = vdASin
    {-# INLINE vASin #-}
    vACos = vdACos
    {-# INLINE vACos #-}
    vATan = vdATan
    {-# INLINE vATan #-}
    vSinh = vdSinh
    {-# INLINE vSinh #-}
    vCosh = vdCosh
    {-# INLINE vCosh #-}
    vTanh = vdTanh
    {-# INLINE vTanh #-}
    vASinh = vdASinh
    {-# INLINE vASinh #-}
    vACosh = vdACosh
    {-# INLINE vACosh #-}
    vATanh = vdATanh
    {-# INLINE vATanh #-}

instance VFloating (Complex Double) where
    vSqrt = vzSqrt
    {-# INLINE vSqrt #-}
    vLog = vzLog
    {-# INLINE vLog #-}

    {- These functions have branch cuts in the wrong places
    vExp = vzExp
    {-# INLINE vExp #-}
    vPow = vzPow
    {-# INLINE vPow #-}
    vSin = vzSin
    {-# INLINE vSin #-}
    vCos = vzCos
    {-# INLINE vCos #-}
    vTan = vzTan
    {-# INLINE vTan #-}
    vASin = vzASin
    {-# INLINE vASin #-}
    vACos = vzACos
    {-# INLINE vACos #-}
    vATan = vzATan
    {-# INLINE vATan #-}
    vSinh = vzSinh
    {-# INLINE vSinh #-}
    vCosh = vzCosh
    {-# INLINE vCosh #-}
    vTanh = vzTanh
    {-# INLINE vTanh #-}
    vASinh = vzASinh
    {-# INLINE vASinh #-}
    vACosh = vzACosh
    {-# INLINE vACosh #-}
    vATanh = vzATanh
    {-# INLINE vATanh #-}
    -}

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE MultiParamTypeClasses, TypeOperators
           , TypeFamilies, UndecidableInstances, CPP
           , FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances  #-}
{-# LANGUAGE DefaultSignatures   #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -Wall #-}
----------------------------------------------------------------------
-- |
-- Module      :   Data.VectorSpace
-- Copyright   :  (c) Conal Elliott and Andy J Gill 2008
-- License     :  BSD3
--
-- Maintainer  :  conal@conal.net, andygill@ku.edu
-- Stability   :  experimental
--
-- Vector spaces
--
-- This version uses associated types instead of fundeps and
-- requires ghc-6.10 or later
----------------------------------------------------------------------

-- NB: I'm attempting to replace fundeps with associated types.  See
-- NewVectorSpace.hs.  Ran into trouble with type equality constraints not
-- getting propagated.  Manuel Ch is looking into it.
--
-- Blocking bug: http://hackage.haskell.org/trac/ghc/ticket/2448

module Data.VectorSpace
  ( module Data.AdditiveGroup
  , VectorSpace(..), (^/), (^*)
  , InnerSpace(..)
  , lerp, linearCombo, magnitudeSq, magnitude, normalized, project
  ) where

import Control.Applicative (liftA2)
import Data.Complex hiding (magnitude)
import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)
import Data.Ratio

import Data.AdditiveGroup
import Data.MemoTrie

import Data.VectorSpace.Generic
import qualified GHC.Generics as Gnrx
import GHC.Generics (Generic, (:*:)(..))

infixr 7 *^

-- | Vector space @v@.
class AdditiveGroup v => VectorSpace v where
  type Scalar v :: *
  type Scalar v = Scalar (VRep v)
  -- | Scale a vector
  (*^) :: Scalar v -> v -> v
  default (*^) :: (Generic v, VectorSpace (VRep v), Scalar (VRep v) ~ Scalar v)
                    => Scalar v -> v -> v
  μ *^ v = Gnrx.to (μ *^ Gnrx.from v :: VRep v)
  {-# INLINE (*^) #-}

infixr 7 <.>

-- | Adds inner (dot) products.
class (VectorSpace v, AdditiveGroup (Scalar v)) => InnerSpace v where
  -- | Inner/dot product
  (<.>) :: v -> v -> Scalar v
  default (<.>) :: (Generic v, InnerSpace (VRep v), Scalar (VRep v) ~ Scalar v)
                    => v -> v -> Scalar v
  v<.>w = (Gnrx.from v :: VRep v) <.> Gnrx.from w
  {-# INLINE (<.>) #-}

infixr 7 ^/
infixl 7 ^*

-- | Vector divided by scalar
(^/) :: (VectorSpace v, s ~ Scalar v, Fractional s) => v -> s -> v
v ^/ s = recip s *^ v
{-# INLINE (^/) #-}

-- | Vector multiplied by scalar
(^*) :: (VectorSpace v, s ~ Scalar v) => v -> s -> v
(^*) = flip (*^)
{-# INLINE (^*) #-}

-- | Linear interpolation between @a@ (when @t==0@) and @b@ (when @t==1@).

-- lerp :: (VectorSpace v, s ~ Scalar v, Num s) => v -> v -> s -> v
lerp :: VectorSpace v => v -> v -> Scalar v -> v
lerp a b t = a ^+^ t *^ (b ^-^ a)
{-# INLINE lerp #-}

-- | Linear combination of vectors
linearCombo :: VectorSpace v => [(v,Scalar v)] -> v
linearCombo ps = sumV [v ^* s | (v,s) <- ps]
{-# INLINE linearCombo #-}

-- | Square of the length of a vector.  Sometimes useful for efficiency.
-- See also 'magnitude'.
magnitudeSq :: (InnerSpace v, s ~ Scalar v) => v -> s
magnitudeSq v = v <.> v
{-# INLINE magnitudeSq #-}

-- | Length of a vector.   See also 'magnitudeSq'.
magnitude :: (InnerSpace v, s ~ Scalar v, Floating s) =>  v -> s
magnitude = sqrt . magnitudeSq
{-# INLINE magnitude #-}

-- | Vector in same direction as given one but with length of one.
-- Divides by zero for the zero vector.
normalized :: (InnerSpace v, s ~ Scalar v, Floating s) =>  v -> v
normalized v = v ^/ magnitude v
{-# INLINE normalized #-}

-- | @project u v@ computes the projection of @v@ onto @u@.
project :: (InnerSpace v, s ~ Scalar v, Fractional s) => v -> v -> v
project u v = ((v <.> u) / magnitudeSq u) *^ u
{-# INLINE project #-}

#define ScalarType(t) \
  instance VectorSpace (t) where \
    { type Scalar (t) = (t) \
    ; (*^) = (*) } ; \
  instance InnerSpace  (t) where (<.>) = (*)

ScalarType(Int)
ScalarType(Integer)
ScalarType(Double)
ScalarType(Float)
ScalarType(CSChar)
ScalarType(CInt)
ScalarType(CShort)
ScalarType(CLong)
ScalarType(CLLong)
ScalarType(CIntMax)
ScalarType(CDouble)
ScalarType(CFloat)

instance Integral a => VectorSpace (Ratio a) where
  type Scalar (Ratio a) = Ratio a
  (*^) = (*)
instance Integral a => InnerSpace  (Ratio a) where (<.>) = (*)

instance (RealFloat v, VectorSpace v) => VectorSpace (Complex v) where
  type Scalar (Complex v) = Scalar v
  s*^(u :+ v) = s*^u :+ s*^v

instance (RealFloat v, InnerSpace v)
     => InnerSpace (Complex v) where
  (u :+ v) <.> (u' :+ v') = (u <.> u') ^+^ (v <.> v')

-- Hm.  The 'RealFloat' constraint is unfortunate here.  It's due to a
-- questionable decision to place 'RealFloat' into the definition of the
-- 'Complex' /type/, rather than in functions and instances as needed.

-- instance (VectorSpace u,VectorSpace v, s ~ Scalar u, s ~ Scalar v)
--          => VectorSpace (u,v) where
--   type Scalar (u,v) = Scalar u
--   s *^ (u,v) = (s*^u,s*^v)

instance ( VectorSpace u, s ~ Scalar u
         , VectorSpace v, s ~ Scalar v )
      => VectorSpace (u,v) where
  type Scalar (u,v) = Scalar u
  s *^ (u,v) = (s*^u,s*^v)

instance ( InnerSpace u, s ~ Scalar u
         , InnerSpace v, s ~ Scalar v )
    => InnerSpace (u,v) where
  (u,v) <.> (u',v') = (u <.> u') ^+^ (v <.> v')

instance ( VectorSpace u, s ~ Scalar u
         , VectorSpace v, s ~ Scalar v
         , VectorSpace w, s ~ Scalar w )
    => VectorSpace (u,v,w) where
  type Scalar (u,v,w) = Scalar u
  s *^ (u,v,w) = (s*^u,s*^v,s*^w)

instance ( InnerSpace u, s ~ Scalar u
         , InnerSpace v, s ~ Scalar v
         , InnerSpace w, s ~ Scalar w )
    => InnerSpace (u,v,w) where
  (u,v,w) <.> (u',v',w') = u<.>u' ^+^ v<.>v' ^+^ w<.>w'

instance ( VectorSpace u, s ~ Scalar u
         , VectorSpace v, s ~ Scalar v
         , VectorSpace w, s ~ Scalar w
         , VectorSpace x, s ~ Scalar x )
    => VectorSpace (u,v,w,x) where
  type Scalar (u,v,w,x) = Scalar u
  s *^ (u,v,w,x) = (s*^u,s*^v,s*^w,s*^x)

instance ( InnerSpace u, s ~ Scalar u
         , InnerSpace v, s ~ Scalar v
         , InnerSpace w, s ~ Scalar w
         , InnerSpace x, s ~ Scalar x )
    => InnerSpace (u,v,w,x) where
  (u,v,w,x) <.> (u',v',w',x') = u<.>u' ^+^ v<.>v' ^+^ w<.>w' ^+^ x<.>x'

-- Standard instances for a functor applied to a vector space.

-- For 'Maybe', Nothing represents 'zeroV'.  Useful for optimization, since
-- we might not be able to test for 'zeroV', e.g., functions and infinite
-- derivative towers.
instance VectorSpace v => VectorSpace (Maybe v) where
  type Scalar (Maybe v) = Scalar v
  (*^) s = fmap (s *^)

-- instance VectorSpace v => VectorSpace (a -> v) where
--   type Scalar (a -> v) = Scalar v
--   (*^) s = fmap (s *^)

-- No 'InnerSpace' instance for @a -> v@.

-- Or the following definition, which is useful for the higher-order
-- shading dsel in Shady (borrowed from Vertigo).

instance VectorSpace v => VectorSpace (a -> v) where
  type Scalar (a -> v) = a -> Scalar v
  (*^) = liftA2 (*^)

instance InnerSpace v => InnerSpace (a -> v) where
  (<.>) = liftA2 (<.>)

instance (HasTrie a, VectorSpace v) => VectorSpace (a :->: v) where
  type Scalar (a :->: v) = Scalar v
  (*^) s = fmap (s *^)

instance InnerSpace a => InnerSpace (Maybe a) where
  -- dotting with zero (vector) yields zero (scalar)
  Nothing <.> _     = zeroV
  _ <.> Nothing     = zeroV
  Just u <.> Just v = u <.> v

--   mu <.> mv = fromMaybe zeroV (liftA2 (<.>) mu mv)

--   (<.>) = (fmap.fmap) (fromMaybe zeroV) (liftA2 (<.>))

instance VectorSpace a => VectorSpace (Gnrx.Rec0 a s) where
  type Scalar (Gnrx.Rec0 a s) = Scalar a
  μ *^ Gnrx.K1 v = Gnrx.K1 $ μ*^v
  {-# INLINE (*^) #-}
instance VectorSpace (f p) => VectorSpace (Gnrx.M1 i c f p) where
  type Scalar (Gnrx.M1 i c f p) = Scalar (f p)
  μ *^ Gnrx.M1 v = Gnrx.M1 $ μ*^v
  {-# INLINE (*^) #-}
instance (VectorSpace (f p), VectorSpace (g p), Scalar (f p) ~ Scalar (g p))
         => VectorSpace ((f :*: g) p) where
  type Scalar ((f:*:g) p) = Scalar (f p)
  μ *^ (x:*:y) = μ*^x :*: μ*^y
  {-# INLINE (*^) #-}

instance InnerSpace a => InnerSpace (Gnrx.Rec0 a s) where
  Gnrx.K1 v <.> Gnrx.K1 w = v<.>w
  {-# INLINE (<.>) #-}
instance InnerSpace (f p) => InnerSpace (Gnrx.M1 i c f p) where
  Gnrx.M1 v <.> Gnrx.M1 w = v<.>w
  {-# INLINE (<.>) #-}
instance ( InnerSpace (f p), InnerSpace (g p)
         , Scalar (f p) ~ Scalar (g p), Num (Scalar (f p)) )
         => InnerSpace ((f :*: g) p) where
  (x:*:y) <.> (ξ:*:υ) = x<.>ξ + y<.>υ
  {-# INLINE (<.>) #-}

max_stars_repo_name: msakai/vector-space

max_issues_repo_name: msakai/vector-space

max_forks_repo_name: msakai/vector-space


text:
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE Trustworthy #-}

{- |
Module      :  Data.Complex.Cyclotomic
Copyright   :  (c) Scott N. Walck 2012-2017
License     :  GPL-3 (see LICENSE)
Maintainer  :  Scott N. Walck <walck@lvc.edu>
Stability   :  experimental

The cyclotomic numbers are a subset of the complex numbers with
the following properties:

     1.  The cyclotomic numbers are represented exactly, enabling exact
     computations and equality comparisons.

     2.  The cyclotomic numbers contain the Gaussian rationals
     (complex numbers of the form 'p' + 'q' 'i' with 'p' and 'q' rational).
     As a consequence, the cyclotomic numbers are a dense subset of the
     complex numbers.

     3.  The cyclotomic numbers contain the square roots of all rational numbers.

     4.  The cyclotomic numbers form a field:  they are closed under addition, subtraction,
     multiplication, and division.

     5.  The cyclotomic numbers contain the sine and cosine of all rational
     multiples of pi.

     6.  The cyclotomic numbers can be thought of as the rational field extended
     with 'n'th roots of unity for arbitrarily large integers 'n'.

     Floating point numbers do not do well with equality comparison:

>(sqrt 2 + sqrt 3)^2 == 5 + 2 * sqrt 6
> -> False

     "Data.Complex.Cyclotomic" represents these numbers exactly, allowing equality comparison:

>(sqrtRat 2 + sqrtRat 3)^2 == 5 + 2 * sqrtRat 6
> -> True

     'Cyclotomic's can be exported as inexact complex numbers using the 'toComplex' function:

>e 6
> -> -e(3)^2
>real $ e 6
> -> 1/2
>imag $ e 6
> -> -1/2*e(12)^7 + 1/2*e(12)^11
>imag (e 6) == sqrtRat 3 / 2
> -> True
>toComplex $ e 6
> -> 0.5000000000000003 :+ 0.8660254037844384

     The algorithms for cyclotomic numbers are adapted from code by
     Martin Schoenert and Thomas Breuer in the GAP project <http://www.gap-system.org/>
     (in particular source files gap4r4\/src\/cyclotom.c and
     gap4r4\/lib\/cyclotom.gi).
-}

module Data.Complex.Cyclotomic
    ( Cyclotomic(..)
    , i
    , e
    , sqrtInteger
    , sqrtRat
    , sinDeg
    , cosDeg
    , sinRev
    , cosRev
    , gaussianRat
    , polarRat
    , polarRatDeg
    , polarRatRev
    , conj
    , real
    , imag
    , isReal
    , isRat
    , isGaussianRat
    , toComplex
    , toReal
    , toRat
    , goldenRatio
    , dft
    , dftInv
    , rootsQuadEq
    , heron
    )
    where

import Data.List
    ( nub
    )
import Data.Ratio
    ( (%)
    , numerator
    , denominator
    )
import Data.Complex
    ( Complex(..)
    , realPart
    )
import qualified Data.Map as M
    ( Map
    , empty
    , singleton
    , lookup
    , keys
    , elems
    , size
    , fromList
    , toList
    , mapKeys
    , filter
    , insertWith
    , delete
    , map
    , unionWith
    , findWithDefault
    , fromListWith
    )
import Math.NumberTheory.Primes.Factorisation
    ( factorise
    )

-- | A cyclotomic number.
data Cyclotomic = Cyclotomic { order  :: Integer
                             , coeffs :: M.Map Integer Rational
                             } deriving (Eq)

-- | @abs@ and @signum@ are partial functions.
--   A cyclotomic number is not guaranteed to have a cyclotomic absolute value.
--   When defined, @signum c@ is the complex number with magnitude 1 that has the same argument as c;
--   @signum c = c / abs c@.
instance Num Cyclotomic where
    (+) = sumCyc
    (*) = prodCyc
    (-) c1 c2 = sumCyc c1 (aInvCyc c2)
    negate = aInvCyc
    abs = absVal
    signum = sigNum
    fromInteger 0 = zeroCyc
    fromInteger n = Cyclotomic 1 (M.singleton 0 (fromIntegral n))

instance Fractional Cyclotomic where
    recip = invCyc
    fromRational 0 = zeroCyc
    fromRational r = Cyclotomic 1 (M.singleton 0 r)

-- | The primitive @n@th root of unity.
--   For example, @'e'(4) = 'i'@ is the primitive 4th root of unity,
--   and 'e'(5) = exp(2*pi*i/5) is the primitive 5th root of unity.
--   In general, 'e' 'n' = exp(2*pi*i/'n').
e :: Integer -> Cyclotomic
e n
    | n < 1      = error "e requires a positive integer"
    | n == 1     = Cyclotomic 1 (M.singleton 0 1)
    | otherwise  = cyclotomic n $ convertToBase n (M.singleton 1 1)

instance Show Cyclotomic where
    show (Cyclotomic n mp)
        | mp == M.empty  = "0"
        | otherwise      = leadingTerm rat n ex ++ followingTerms n xs
        where ((ex,rat):xs) = M.toList mp

showBaseExp :: Integer -> Integer -> String
showBaseExp n 1  = "e(" ++ show n ++ ")"
showBaseExp n ex = "e(" ++ show n ++ ")^" ++ show ex

leadingTerm :: Rational -> Integer -> Integer -> String
leadingTerm r _ 0 = showRat r
leadingTerm r n ex
    | r == 1     = t
    | r == (-1)  = '-' : t
    | r > 0      = showRat r ++ "*" ++ t
    | r < 0      = "-" ++ showRat (abs r) ++ "*" ++ t
    | otherwise  = ""
    where t = showBaseExp n ex

followingTerms :: Integer -> [(Integer,Rational)] -> String
followingTerms _ [] = ""
followingTerms n ((ex,rat):xs) = followingTerm rat n ex ++ followingTerms n xs

followingTerm :: Rational -> Integer -> Integer -> String
followingTerm r n ex
    | r == 1     = " + " ++ t
    | r == (-1)  = " - " ++ t
    | r > 0      = " + " ++ showRat r ++ "*" ++ t
    | r < 0      = " - " ++ showRat (abs r) ++ "*" ++ t
    | otherwise  = ""
    where t = showBaseExp n ex

showRat :: Rational -> String
showRat r
    | d == 1     = show n
    | otherwise  = show n ++ "/" ++ show d
    where
      n = numerator r
      d = denominator r

-- GAP function EB from gap4r4/lib/cyclotom.gi
eb :: Integer -> Cyclotomic
eb n
    | n < 1           = error "eb needs a positive integer"
    | n `mod` 2 /= 1  = error "eb needs an odd integer"
    | n == 1          = zeroCyc
    | otherwise       = let en = e n
                        in sum [en^(k*k `mod` n) | k <- [1..(n-1) `div` 2]]

sqrt2 :: Cyclotomic
sqrt2 = e 8 - e 8 ^ (3 :: Int)

-- | The square root of an 'Integer'.
sqrtInteger :: Integer -> Cyclotomic
sqrtInteger n
    | n == 0     = zeroCyc
    | n < 0      = i * sqrtPositiveInteger (-n)
    | otherwise  = sqrtPositiveInteger n

sqrtPositiveInteger :: Integer -> Cyclotomic
sqrtPositiveInteger n
    | n < 1      = error "sqrtPositiveInteger needs a positive integer"
    | otherwise  = let factors = factorise n
                       factor = product [p^(m `div` 2) | (p,m) <- factors]
                       nn     = product [p^(m `mod` 2) | (p,m) <- factors]
                   in case nn `mod` 4 of
                        1 -> fromInteger factor * (2 * eb nn + 1)
                        2 -> fromInteger factor * sqrt2 * sqrtPositiveInteger (nn `div` 2)
                        3 -> fromInteger factor * (-i) * (2 * eb nn + 1)
                        _ -> fromInteger factor * 2 * sqrtPositiveInteger (nn `div` 4)

-- | The square root of a 'Rational' number.
sqrtRat :: Rational -> Cyclotomic
sqrtRat r = prodRatCyc (1 % fromInteger den) (sqrtInteger (numerator r * den))
    where
      den = denominator r

-- | The square root of -1.
i :: Cyclotomic
i = e 4

-- | Make a Gaussian rational; @gaussianRat p q@ is the same as @p + q * i@.
gaussianRat :: Rational -> Rational -> Cyclotomic
gaussianRat p q = fromRational p + fromRational q * i

-- | A complex number in polar form, with rational magnitude @r@ and rational angle @s@
--   of the form @r * exp(2*pi*i*s)@; @polarRat r s@ is the same as @r * e q ^ p@,
--   where @s = p/q@.  This function is the same as 'polarRatRev'.
polarRat :: Rational    -- ^ magnitude
         -> Rational    -- ^ angle, in revolutions
         -> Cyclotomic  -- ^ cyclotomic number
polarRat r s
    = let p = numerator s
          q = denominator s
      in case p >= 0 of
           True  -> fromRational r * e q ^ p
           False -> conj $ fromRational r * e q ^ (-p)

-- | A complex number in polar form, with rational magnitude and rational angle
--   in degrees.
polarRatDeg :: Rational    -- ^ magnitude
            -> Rational    -- ^ angle, in degrees
            -> Cyclotomic  -- ^ cyclotomic number
polarRatDeg r deg
    = let s = deg / 360
          p = numerator s
          q = denominator s
      in case p >= 0 of
           True  -> fromRational r * e q ^ p
           False -> conj $ fromRational r * e q ^ (-p)

-- | A complex number in polar form, with rational magnitude and rational angle
--   in revolutions.
polarRatRev :: Rational    -- ^ magnitude
            -> Rational    -- ^ angle, in revolutions
            -> Cyclotomic  -- ^ cyclotomic number
polarRatRev r s
    = let p = numerator s
          q = denominator s
      in case p >= 0 of
           True  -> fromRational r * e q ^ p
           False -> conj $ fromRational r * e q ^ (-p)

-- | Complex conjugate.
conj :: Cyclotomic -> Cyclotomic
conj (Cyclotomic n mp)
    = mkCyclotomic n (M.mapKeys (\k -> (n-k) `mod` n) mp)

-- | Real part of the cyclotomic number.
real :: Cyclotomic -> Cyclotomic
real z = (z + conj z) / 2

-- | Imaginary part of the cyclotomic number.
imag :: Cyclotomic -> Cyclotomic
imag z = (z - conj z) / (2*i)

absVal :: Cyclotomic -> Cyclotomic
absVal c = let modsq = c * conj c
           in case toRat modsq of
                Just msq -> sqrtRat msq
                Nothing  -> error "abs not available for this number"

sigNum :: Cyclotomic -> Cyclotomic
sigNum 0 = zeroCyc
sigNum c = let modsq = c * conj c
           in if isRat modsq then c / absVal c
              else error "signum not available for this number"

convertToBase :: Integer -> M.Map Integer Rational -> M.Map Integer Rational
convertToBase n mp = foldr (\(p,r) m -> replace n p r m) mp (extraneousPowers n)

removeZeros :: M.Map Integer Rational -> M.Map Integer Rational
removeZeros = M.filter (/= 0)

-- Corresponds to GAP implementation.
-- Expects that convertToBase has already been done.
cyclotomic :: Integer -> M.Map Integer Rational -> Cyclotomic
cyclotomic ord = tryReduce . tryRational . gcdReduce . Cyclotomic ord

mkCyclotomic :: Integer -> M.Map Integer Rational -> Cyclotomic
mkCyclotomic ord = cyclotomic ord . removeZeros . convertToBase ord

-- | Step 1 of cyclotomic is gcd reduction.
gcdReduce :: Cyclotomic -> Cyclotomic
gcdReduce cyc@(Cyclotomic n mp) = case gcdCyc cyc of
                                    1 -> cyc
                                    d -> Cyclotomic (n `div` d) (M.mapKeys (`div` d) mp)

gcdCyc :: Cyclotomic -> Integer
gcdCyc (Cyclotomic n mp) = gcdList (n:M.keys mp)

-- | Step 2 of cyclotomic is reduction to a rational if possible.
tryRational :: Cyclotomic -> Cyclotomic
tryRational c
    | lenCyc c == fromIntegral phi && sqfree
        = case equalCoefficients c of
            Nothing -> c
            Just r  -> fromRational $ (-1)^(nrp `mod` 2)*r
    | otherwise
        = c
    where
      (phi,nrp,sqfree) = phiNrpSqfree (order c)

-- | Compute phi(n), the number of prime factors, and test if n is square-free.
--   We do these all together for efficiency, so we only call factorise once.
phiNrpSqfree :: Integer -> (Integer,Int,Bool)
phiNrpSqfree n = (phi,nrp,sqfree)
    where
      factors = factorise n
      phi = foldr (\p n' -> n' `div` p * (p-1)) n [p | (p,_) <- factors]
      nrp = length factors
      sqfree = all (<=1) [m | (_,m) <- factors]

equalCoefficients :: Cyclotomic -> Maybe Rational
equalCoefficients (Cyclotomic _ mp)
    = case ts of
        []    -> Nothing
        (x:_) -> if equal ts then Just x else Nothing
      where
        ts = M.elems mp

lenCyc :: Cyclotomic -> Int
lenCyc (Cyclotomic _ mp) = M.size $ removeZeros mp

-- | Step 3 of cyclotomic is base reduction
tryReduce :: Cyclotomic -> Cyclotomic
tryReduce c
    = foldr reduceByPrime c squareFreeOddFactors
      where
        squareFreeOddFactors = [p | (p,m) <- factorise (order c), p > 2, m <= 1]

reduceByPrime :: Integer -> Cyclotomic -> Cyclotomic
reduceByPrime p c@(Cyclotomic n _)
    = case mapM (\r -> equalReplacements p r c) [0,p..n-p] of
        Just cfs -> Cyclotomic (n `div` p) $ removeZeros $ M.fromList $ zip [0..(n `div` p)-1] (map negate cfs)
        Nothing  -> c

equalReplacements :: Integer -> Integer -> Cyclotomic -> Maybe Rational
equalReplacements p r (Cyclotomic n mp)
    =  case [M.findWithDefault 0 k mp | k <- replacements n p r] of
         [] -> error "equalReplacements generated empty list"
         (x:xs) | equal (x:xs) -> Just x
         _ -> Nothing

replacements :: Integer -> Integer -> Integer -> [Integer]
replacements n p r = takeWhile (>= 0) [r-s,r-2*s..] ++ takeWhile (< n) [r+s,r+2*s..]
    where s = n `div` p

replace :: Integer -> Integer -> Integer -> M.Map Integer Rational -> M.Map Integer Rational
replace n p r mp = case M.lookup r mp of
                     Nothing  -> mp
                     Just rat -> foldr (\k m -> M.insertWith (+) k (-rat) m) (M.delete r mp) (replacements n p r)

includeMods :: Integer -> Integer -> Integer -> [Integer]
includeMods n q start = [start] ++ takeWhile (>= 0) [start-q,start-2*q..] ++ takeWhile (< n) [start+q,start+2*q..]

removeExps :: Integer -> Integer -> Integer -> [Integer]
removeExps n 2 q = concatMap (includeMods n q) $ map ((n `div` q) *) [q `div` 2..q-1]
removeExps n p q = concatMap (includeMods n q) $ map ((n `div` q) *) [-m..m]
    where m = (q `div` p - 1) `div` 2

pqPairs :: Integer -> [(Integer,Integer)]
pqPairs n = map (\(p,k) -> (p,p^k)) (factorise n)

extraneousPowers :: Integer -> [(Integer,Integer)]
extraneousPowers n
    | n < 1      = error "extraneousPowers needs a postive integer"
    | otherwise  = nub $ concat [[(p,r) | r <- removeExps n p q] | (p,q) <- pqPairs n]

-- | Sum of two cyclotomic numbers.
sumCyc :: Cyclotomic -> Cyclotomic -> Cyclotomic
sumCyc (Cyclotomic o1 map1) (Cyclotomic o2 map2)
    = let ord = lcm o1 o2
          m1 = ord `div` o1
          m2 = ord `div` o2
          map1' = M.mapKeys (m1*) map1
          map2' = M.mapKeys (m2*) map2
      in mkCyclotomic ord $ M.unionWith (+) map1' map2'

-- | Product of two cyclotomic numbers.
prodCyc :: Cyclotomic -> Cyclotomic -> Cyclotomic
prodCyc (Cyclotomic o1 map1) (Cyclotomic o2 map2)
    = let ord = lcm o1 o2
          m1 = ord `div` o1
          m2 = ord `div` o2
      in mkCyclotomic ord $ M.fromListWith (+)
             [((m1*e1+m2*e2) `mod` ord,c1*c2) | (e1,c1) <- M.toList map1, (e2,c2) <- M.toList map2]

-- | Product of a rational number and a cyclotomic number.
prodRatCyc :: Rational -> Cyclotomic -> Cyclotomic
prodRatCyc 0 _                   = zeroCyc
prodRatCyc r (Cyclotomic ord mp) = Cyclotomic ord $ M.map (r*) mp

-- | Additive identity.
zeroCyc :: Cyclotomic
zeroCyc = Cyclotomic 1 M.empty

-- | Additive inverse.
aInvCyc :: Cyclotomic -> Cyclotomic
aInvCyc = prodRatCyc (-1)

multiplyExponents :: Integer -> Cyclotomic -> Cyclotomic
multiplyExponents j (Cyclotomic n mp)
    | gcd j n /= 1  = error "multiplyExponents needs gcd == 1"
    | otherwise     = mkCyclotomic n (M.mapKeys (\k -> j*k `mod` n) mp)

productOfGaloisConjugates :: Cyclotomic -> Cyclotomic
productOfGaloisConjugates c
    = product [multiplyExponents j c | j <- [2..ord], gcd j ord == 1]
      where
        ord = order c

-- | Multiplicative inverse.
invCyc :: Cyclotomic -> Cyclotomic
invCyc z = case toRat (z * prod) of
             Just r  -> prodRatCyc (1 / r) prod
             Nothing -> error "invCyc:  product of Galois conjugates not rational; this is a bug, please inform package maintainer"
    where
      prod = productOfGaloisConjugates z

-- | Is the cyclotomic a real number?
isReal :: Cyclotomic -> Bool
isReal c = c == conj c

-- | Is the cyclotomic a rational?
isRat :: Cyclotomic -> Bool
isRat (Cyclotomic 1 _) = True
isRat _                = False

-- | Is the cyclotomic a Gaussian rational?
isGaussianRat :: Cyclotomic -> Bool
isGaussianRat c = isRat (real c) && isRat (imag c)

-- | Original implementation of toComplex
origToComplex :: RealFloat a => Cyclotomic -> Complex a
origToComplex c = sum [fromRational r * en^p | (p,r) <- M.toList (coeffs c)]
    where en = exp (0 :+ 2*pi/n)
          n = fromIntegral (order c)

-- | Export as an inexact complex number.
-- This function avoids some error that origToComplex introduces.
toComplex :: RealFloat a => Cyclotomic -> Complex a
toComplex x = re :+ im
  where
    Just re = fromRational . toRational <$> (toReal (real x) :: Maybe Double)
    Just im = fromRational . toRational <$> (toReal (imag x) :: Maybe Double)

-- | Export as an inexact real number if possible.
toReal :: RealFloat a => Cyclotomic -> Maybe a
toReal c
    | isReal c   = Just $ realPart (origToComplex c)
    | otherwise  = Nothing

-- | Return an exact rational number if possible.
toRat :: Cyclotomic -> Maybe Rational
toRat (Cyclotomic 1 mp)
    | mp == M.empty  = Just 0
    | otherwise      = M.lookup 0 mp
toRat _ = Nothing

-- | Sine function with argument in degrees.
sinDeg :: Rational -> Cyclotomic
sinDeg d = let n = d / 360
               nm = abs (numerator n)
               dn = denominator n
               a = e dn^nm
           in fromRational(signum d) * (a - conj a) / (2*i)

-- | Cosine function with argument in degrees.
cosDeg :: Rational -> Cyclotomic
cosDeg d = let n = d / 360
               nm = abs (numerator n)
               dn = denominator n
               a = e dn^nm
           in (a + conj a) / 2

-- | Sine function with argument in revolutions.
sinRev :: Rational -> Cyclotomic
sinRev n = let nm = abs (numerator n)
               dn = denominator n
               a = e dn^nm
           in fromRational(signum n) * (a - conj a) / (2*i)

-- | Cosine function with argument in revolutions.
cosRev :: Rational -> Cyclotomic
cosRev n = let nm = abs (numerator n)
               dn = denominator n
               a = e dn^nm
           in (a + conj a) / 2

gcdList :: [Integer] -> Integer
gcdList [] = error "gcdList called on empty list"
gcdList (n:ns) = foldr gcd n ns

equal :: Eq a => [a] -> Bool
equal [] = True
equal [_] = True
equal (x:y:ys) = x == y && equal (y:ys)

-- | The golden ratio, @(1 + √5)/2@.
goldenRatio :: Cyclotomic
goldenRatio = (1 + sqrtRat 5) / 2

-- | Discrete Fourier transform,
--   @X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i 2 \pi \frac{k}{N} n}@.
dft :: [Cyclotomic] -> [Cyclotomic]
dft cs = [sum $ zipWith (*) [conj (e m^(k*n)) | n <- [0..]] cs | k <- [0..m-1]]
          where m = fromIntegral $ length cs

-- | Inverse discrete Fourier transform,
--   @x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k \cdot e^{i 2 \pi \frac{k}{N} n}@.
dftInv :: [Cyclotomic] -> [Cyclotomic]
dftInv cs = [minv * sum (zipWith (*) [e m^(k*n) | n <- [0..]] cs) | k <- [0..m-1]]
          where m = fromIntegral $ length cs
                minv = fromRational (1 % m)

-- | Solutions to the quadratic equation
--   a x^2 + b x + c = 0.
--   Returns 'Nothing' if a == 0.
rootsQuadEq :: Rational  -- ^ a
            -> Rational  -- ^ b
            -> Rational  -- ^ c
            -> Maybe (Cyclotomic,Cyclotomic)  -- ^ roots
rootsQuadEq a b c
    | a == 0     = Nothing
    | otherwise  = Just ((-bb + sqrtDisc)/(2*aa),(-bb - sqrtDisc)/(2*aa))
    where
      aa = fromRational a
      bb = fromRational b
      sqrtDisc = sqrtRat (b*b - 4*a*c)

-- | Heron's formula for the area of a triangle with
--   side lengths a, b, c.
heron :: Rational    -- ^ a
      -> Rational    -- ^ b
      -> Rational    -- ^ c
      -> Cyclotomic  -- ^ area of triangle
heron a b c
    = sqrtRat (s * (s-a) * (s-b) * (s-c))
      where
        s = (a + b + c) / 2

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
{-# LANGUAGE CPP                       #-}
{-# LANGUAGE ConstraintKinds           #-}
{-# LANGUAGE FlexibleContexts          #-}

module TSA

where

import Control.Monad.IO.Class
import Data.Attoparsec.ByteString.Char8
import Data.CSV.Conduit
import Data.Maybe
#if MIN_VERSION_base(4,11,0)
#else
import Data.Monoid
#endif

import qualified Data.Vector.Unboxed as V
import qualified Data.Vector as VB
import qualified Statistics.Regression as S
import qualified Statistics.Sample as S

#ifdef WITH_IHASKELL
import Control.Lens
import Graphics.Rendering.Chart.Easy as Chart

import qualified Statistics.Autocorrelation as S
import qualified Statistics.Sample.Histogram as S
#endif

type Series e = V.Vector e

loadSeries :: MonadIO m
           => FilePath
           -- ^ CSV file name.
           -> Int
           -- ^ Index of column to use (1-based).
           -> m (Series Double)
loadSeries fname col = do
  r <- readCSVFile defCSVSettings fname
  return $ V.convert $ flip VB.mapMaybe r $
    \row -> (either (const Nothing) Just . parseOnly double) =<< row VB.!? col

diff :: (Num e, V.Unbox e)
     => Int
     -- ^ Differencing period.
     -> Series e
     -> Series e
diff n series = V.zipWith (-) (V.drop n series) series

-- | Model is a forecasting function.
--
-- If the argument is Nothing, predict existing values.
type Model e = Maybe Int -> Series e -> Series e
{-@ type Model' e = Maybe Nat -> Series e -> Series e @-}

-- | Predict N next points for a time series.
{-@ forecast :: Nat -> Model' e -> Series e -> Series e @-}
forecast :: Int -> Model e -> Series e -> Series e
forecast periods f = f $ Just periods

{-@ extendBy :: Nat -> Model' e -> Series e -> Series e @-}
extendBy :: V.Unbox e => Int -> Model e -> Series e -> Series e
extendBy periods f s = s <> forecast periods f s

{-@ mean :: Model' Double @-}
mean :: Model Double
mean n s = V.replicate n' $ S.mean s
  where
    n' = fromMaybe (V.length s) n

{-@ linearTrend :: Model' Double @-}
linearTrend :: Model Double
linearTrend n s = V.iterateN n' (+ slope) start
  where
    n' = fromMaybe len n
    len = V.length s
    (coeffs, _) = S.olsRegress [V.generate len fromIntegral] s
    (slope, bias) = (coeffs V.! 0, coeffs V.! 1)
    start = fromIntegral (maybe 0 (const len) n) * slope + bias

residuals :: (Num e, V.Unbox e) => Model e -> Series e -> Series e
residuals f s = V.zipWith (-) s (f Nothing s)

#ifdef WITH_IHASKELL
plotResiduals :: Model Double -> Series Double -> Renderable ()
plotResiduals f s = toRenderable $ do
  layout_title .= "Residuals"
  layout_legend .= Nothing
  Chart.plot $ plotBars <$> do
    b <- bars (map show $ V.toList sizes) (map (\(x, y) -> (x, [y :: Int])) barData)
    return $ b & plot_bars_spacing .~ BarsFixGap 0 5
  where
    barData = V.toList $ V.zip bins sizes
    -- TODO Customize bin count?
    (bins, sizes) = S.histogram 20 $ residuals f s

plotACF :: Int -> Series Double -> Renderable ()
plotACF maxLag s = toRenderable $ do
  layout_title .= "Autocorrelation"
  layout_legend .= Nothing
  layout_x_axis . laxis_title .= "Lag"
  -- Do not skip X values
  layout_x_axis . laxis_override .= (axis_labels .~ [[ (lag, show lag) | lag <- lags ]])
  Chart.plot $ plotBars <$>
    bars (map show lags) (map (\lag -> (lag, [acf V.! lag])) lags)
  where
    (acf, _, _) = S.autocorrelation s
    lags = [1 .. maxLag]

plotPredictions :: (PlotValue e, V.Unbox e)
                => Int
                -- ^ Forecast that many future points.
                -> Model e
                -> Series e
                -> Renderable ()
plotPredictions periods f s = toRenderable $ do
  Chart.setColors [opaque blue, red `withOpacity` 50, opaque red]
  Chart.plot $ line "Series" [zip [(0::Int)..] (V.toList s)]
  mapM_ (plotWithRadius 5)
    [ points "Predictions" $ zip [(0::Int)..] (V.toList $ f Nothing s)
    , points "Forecast" $ zip [V.length s..] (V.toList $ forecast periods f s)
    ]
  where
    plotWithRadius radius p = Chart.plot $ do
      p' <- p
      return (p' & plot_points_style . point_radius .~ radius)
#endif

max_stars_repo_name: dzhus/tsa

max_issues_repo_name: dzhus/tsa

max_forks_repo_name: dzhus/tsa


text:
-- | See Hoffman, Gelman (2011) The No U-Turn Sampler: Adaptively Setting Path
--   Lengths in Hamiltonian Monte Carlo.

module Numeric.MCMC.NUTS where

import Control.Monad
import Control.Monad.Loops
import Control.Monad.Primitive
import System.Random.MWC
import System.Random.MWC.Distributions hiding (gamma)
import Statistics.Distribution.Normal

type Parameters = [Double]
type Density    = Parameters -> Double
type Gradient   = Parameters -> Parameters
type Particle   = (Parameters, Parameters)
type StateSpecs = ([Double], [Double], [Double], [Double], [Double], Int, Int)

newtype StateTree = StateTree {
    getStateTree :: StateSpecs
  }

instance Show StateTree where
  show (StateTree (tm, rm, tp, rp, t', n, s)) =
       "\n" ++ "tm: " ++ show tm
    ++ "\n" ++ "tp: " ++ show tp
    ++ "\n" ++ "t': " ++ show t'
    ++ "\n" ++ "n : " ++ show n
    ++ "\n" ++ "s : " ++ show s

data DualAveragingParameters = DualAveragingParameters {
    mAdapt :: Int
  , delta  :: Double
  , mu     :: Double
  , gamma  :: Double
  , tau0   :: Double
  , kappa  :: Double
  } deriving Show

-- | The NUTS sampler.
nuts
  :: PrimMonad m
  => Density
  -> Gradient
  -> Int
  -> Double
  -> Parameters
  -> Gen (PrimState m)
  -> m [Parameters]
nuts lTarget glTarget n e t g = go t 0 []
  where go position j acc
          | j >= n    = return acc
          | otherwise = do
              nextPosition <- nutsKernel lTarget glTarget e position g
              go nextPosition (succ j) (nextPosition : acc)

-- | The NUTS sampler with dual averaging.
nutsDualAveraging
  :: PrimMonad m
  => Density
  -> Gradient
  -> Int
  -> Int
  -> Parameters
  -> Gen (PrimState m)
  -> m [Parameters]
nutsDualAveraging lTarget glTarget n nAdapt t g = do
    e0 <- findReasonableEpsilon lTarget glTarget t g
    let daParams = basicDualAveragingParameters e0 nAdapt
    chain <- unfoldrM (kernel daParams) (1, e0, 1, 0, t)
    return $ drop nAdapt chain
  where
    kernel params (m, e, eAvg, h, t0) = do
      (eNext, eAvgNext, hNext, tNext) <-
        nutsKernelDualAvg lTarget glTarget e eAvg h m params t0 g
      return $ if   m > n + nAdapt
               then Nothing
               else Just (t0, (succ m, eNext, eAvgNext, hNext, tNext))

-- | Default DA parameters, given a base step size and burn in period.
basicDualAveragingParameters :: Double -> Int -> DualAveragingParameters
basicDualAveragingParameters step burnInPeriod = DualAveragingParameters {
    mu     = log (10 * step)
  , delta  = 0.5
  , mAdapt = burnInPeriod
  , gamma  = 0.05
  , tau0   = 10
  , kappa  = 0.75
  }

-- | A single iteration of dual-averaging NUTS.
nutsKernelDualAvg
  :: PrimMonad m
  => Density
  -> Gradient
  -> Double
  -> Double
  -> Double
  -> Int
  -> DualAveragingParameters
  -> [Double]
  -> Gen (PrimState m)
  -> m (Double, Double, Double, Parameters)
nutsKernelDualAvg lTarget glTarget e eAvg h m daParams t g = do
  r0 <- replicateM (length t) (normal 0 1 g)
  z0 <- exponential 1 g
  let logu = log (auxilliaryTarget lTarget t r0) - z0

  let go (tn, tp, rn, rp, tm, j, n, s, a, na) g
        | s == 1 = do
            vj <- symmetricCategorical [-1, 1] g
            z  <- uniform g

            (tnn, rnn, tpp, rpp, t1, n1, s1, a1, na1) <-
              if   vj == -1
              then do
                (tnn', rnn', _, _, t1', n1', s1', a1', na1') <-
                  buildTreeDualAvg lTarget glTarget g tn rn logu vj j e t r0
                return (tnn', rnn', tp, rp, t1', n1', s1', a1', na1')
              else do
                (_, _, tpp', rpp', t1', n1', s1', a1', na1') <-
                  buildTreeDualAvg lTarget glTarget g tp rp logu vj j e t r0
                return (tn, rn, tpp', rpp', t1', n1', s1', a1', na1')

            let accept = s1 == 1 && (min 1 (fi n1 / fi n :: Double)) > z

                n2 = n + n1
                s2 = s1 * stopCriterion tnn tpp rnn rpp
                j1 = succ j
                t2 | accept    = t1
                   | otherwise = tm

            go (tnn, tpp, rnn, rpp, t2, j1, n2, s2, a1, na1) g

        | otherwise = return (tm, a, na)

  (nextPosition, alpha, nalpha) <- go (t, t, r0, r0, t, 0, 1, 1, 0, 0) g

  let (hNext, eNext, eAvgNext) =
          if   m <= mAdapt daParams
          then (hm, exp logEm, exp logEbarM)
          else (h, eAvg, eAvg)
        where
          eta = 1 / (fromIntegral m + tau0 daParams)
          hm  = (1 - eta) * h
              + eta * (delta daParams - alpha / fromIntegral nalpha)

          zeta = fromIntegral m ** (- (kappa daParams))

          logEm    = mu daParams - sqrt (fromIntegral m) / gamma daParams * hm
          logEbarM = (1 - zeta) * log eAvg + zeta * logEm

  return (eNext, eAvgNext, hNext, nextPosition)

-- | A single iteration of NUTS.
nutsKernel
  :: PrimMonad m
  => Density
  -> Gradient
  -> Double
  -> Parameters
  -> Gen (PrimState m)
  -> m Parameters
nutsKernel lTarget glTarget e t g = do
  r0   <- replicateM (length t) (normal 0 1 g)
  z0   <- exponential 1 g
  let logu = log (auxilliaryTarget lTarget t r0) - z0

  let go (tn, tp, rn, rp, tm, j, n, s) g
        | s == 1 = do
            vj <- symmetricCategorical [-1, 1] g
            z  <- uniform g

            (tnn, rnn, tpp, rpp, t1, n1, s1) <-
              if   vj == -1
              then do
                (tnn', rnn', _, _, t1', n1', s1') <-
                  buildTree lTarget glTarget g tn rn logu vj j e
                return (tnn', rnn', tp, rp, t1', n1', s1')
              else do
                (_, _, tpp', rpp', t1', n1', s1') <-
                  buildTree lTarget glTarget g tp rp logu vj j e
                return (tn, rn, tpp', rpp', t1', n1', s1')

            let accept = s1 == 1 && (min 1 (fi n1 / fi n :: Double)) > z

                n2 = n + n1
                s2 = s1 * stopCriterion tnn tpp rnn rpp
                j1 = succ j
                t2 | accept    = t1
                   | otherwise = tm

            go (tnn, tpp, rnn, rpp, t2, j1, n2, s2) g

        | otherwise = return tm

  go (t, t, r0, r0, t, 0, 1, 1) g

-- | Build the 'tree' of candidate states.
buildTree
  :: PrimMonad m
  => Density
  -> Gradient
  -> Gen (PrimState m)
  -> Parameters
  -> Parameters
  -> Double
  -> Double
  -> Int
  -> Double
  -> m StateSpecs
buildTree lTarget glTarget g t r logu v 0 e = do
  let (t0, r0) = leapfrog glTarget (t, r) (v * e)
      joint    = log $ auxilliaryTarget lTarget t0 r0
      n        = indicate (logu < joint)
      s        = indicate (logu - 1000 < joint)
  return (t0, r0, t0, r0, t0, n, s)

buildTree lTarget glTarget g t r logu v j e = do
  z <- uniform g
  (tn, rn, tp, rp, t0, n0, s0) <-
    buildTree lTarget glTarget g t r logu v (pred j) e

  if   s0 == 1
  then do
    (tnn, rnn, tpp, rpp, t1, n1, s1) <-
      if   v == -1
      then do
        (tnn', rnn', _, _, t1', n1', s1') <-
          buildTree lTarget glTarget g tn rn logu v (pred j) e
        return (tnn', rnn', tp, rp, t1', n1', s1')
      else do
        (_, _, tpp', rpp', t1', n1', s1') <-
          buildTree lTarget glTarget g tp rp logu v (pred j) e
        return (tn, rn, tpp', rpp', t1', n1', s1')

    let accept = (fi n1 / max (fi (n0 + n1)) 1) > (z :: Double)
        n2     = n0 + n1
        s2     = s0 * s1 * stopCriterion tnn tpp rnn rpp
        t2     | accept    = t1
               | otherwise = t0

    return (tnn, rnn, tpp, rpp, t2, n2, s2)
  else return (tn, rn, tp, rp, t0, n0, s0)

-- | Determine whether or not to stop doubling the tree of candidate states.
stopCriterion :: (Integral a, Num b, Ord b) => [b] -> [b] -> [b] -> [b] -> a
stopCriterion tn tp rn rp =
      indicate (positionDifference `innerProduct` rn >= 0)
    * indicate (positionDifference `innerProduct` rp >= 0)
  where
    positionDifference = tp .- tn

-- | Build the tree of candidate states under dual averaging.
buildTreeDualAvg
  :: PrimMonad m
  => Density
  -> Gradient
  -> Gen (PrimState m)
  -> Parameters
  -> Parameters
  -> Double
  -> Double
  -> Int
  -> Double
  -> Parameters
  -> Parameters
  -> m ([Double], [Double], [Double], [Double], [Double], Int, Int, Double, Int)
buildTreeDualAvg lTarget glTarget g t r logu v 0 e t0 r0 = do
  let (t1, r1) = leapfrog glTarget (t, r) (v * e)
      joint    = log $ auxilliaryTarget lTarget t1 r1
      n        = indicate (logu < joint)
      s        = indicate (logu - 1000 <  joint)
      a        = min 1 (acceptanceRatio lTarget t0 t1 r0 r1)
  return (t1, r1, t1, r1, t1, n, s, a, 1)

buildTreeDualAvg lTarget glTarget g t r logu v j e t0 r0 = do
  z <- uniform g
  (tn, rn, tp, rp, t1, n1, s1, a1, na1) <-
    buildTreeDualAvg lTarget glTarget g t r logu v (pred j) e t0 r0

  if   s1 == 1
  then do
    (tnn, rnn, tpp, rpp, t2, n2, s2, a2, na2) <-
      if   v == -1
      then do
        (tnn', rnn', _, _, t1', n1', s1', a1', na1') <-
          buildTreeDualAvg lTarget glTarget g tn rn logu v (pred j) e t0 r0
        return (tnn', rnn', tp, rp, t1', n1', s1', a1', na1')
      else do
        (_, _, tpp', rpp', t1', n1', s1', a1', na1') <-
          buildTreeDualAvg lTarget glTarget g tp rp logu v (pred j) e t0 r0
        return (tn, rn, tpp', rpp', t1', n1', s1', a1', na1')

    let p      = fi n2 / max (fi (n1 + n2)) 1
        accept = p > (z :: Double)
        n3     = n1 + n2
        a3     = a1 + a2
        na3    = na1 + na2
        s3     = s1 * s2 * stopCriterion tnn tpp rnn rpp

        t3  | accept    = t2
            | otherwise = t1

    return (tnn, rnn, tpp, rpp, t3, n3, s3, a3, na3)
  else return (tn, rn, tp, rp, t1, n1, s1, a1, na1)

-- | Heuristic for initializing step size.
findReasonableEpsilon
  :: PrimMonad m
  => Density
  -> Gradient
  -> Parameters
  -> Gen (PrimState m)
  -> m Double
findReasonableEpsilon lTarget glTarget t0 g = do
  r0 <- replicateM (length t0) (normal 0 1 g)
  let (t1, r1) = leapfrog glTarget (t0, r0) 1.0
      a        = 2 * indicate (acceptanceRatio lTarget t0 t1 r0 r1 > 0.5) - 1

      go j e t r
        | j <= 0 = e -- no need to shrink this excessively
        | (acceptanceRatio lTarget t0 t r0 r) ^^ a > 2 ^^ (-a) =
            let (tn, rn) = leapfrog glTarget (t, r) e
            in  go (pred j) (2 ^^ a * e) tn rn
        | otherwise = e

  return $ go 10 1.0 t1 r1

-- | Simulate a single step of Hamiltonian dynamics.
leapfrog :: Gradient -> Particle -> Double -> Particle
leapfrog glTarget (t, r) e = (tf, rf)
  where
    rm = adjustMomentum glTarget e t r
    tf = adjustPosition e rm t
    rf = adjustMomentum glTarget e tf rm

-- | Adjust momentum.
adjustMomentum :: Fractional c => (t -> [c]) -> c -> t -> [c] -> [c]
adjustMomentum glTarget e t r = r .+ ((e / 2) .* glTarget t)

-- | Adjust position.
adjustPosition :: Num c => c -> [c] -> [c] -> [c]
adjustPosition e r t = t .+ (e .* r)

-- | The MH acceptance ratio for a given proposal.
acceptanceRatio :: Floating a => (t -> a) -> t -> t -> [a] -> [a] -> a
acceptanceRatio lTarget t0 t1 r0 r1 = auxilliaryTarget lTarget t1 r1
                                    / auxilliaryTarget lTarget t0 r0

-- | The negative potential.
auxilliaryTarget :: Floating a => (t -> a) -> t -> [a] -> a
auxilliaryTarget lTarget t r = exp (lTarget t - 0.5 * innerProduct r r)

-- | Simple inner product.
innerProduct :: Num a => [a] -> [a] -> a
innerProduct xs ys = sum $ zipWith (*) xs ys

-- | Vectorized multiplication.
(.*) :: Num b => b -> [b] -> [b]
z .* xs = map (* z) xs

-- | Vectorized subtraction.
(.-) :: Num a => [a] -> [a] -> [a]
xs .- ys = zipWith (-) xs ys

-- | Vectorized addition.
(.+) :: Num a => [a] -> [a] -> [a]
xs .+ ys = zipWith (+) xs ys

-- | Indicator function.
indicate :: Integral a => Bool -> a
indicate True  = 1
indicate False = 0

-- | A symmetric categorical (discrete uniform) distribution.
symmetricCategorical :: PrimMonad m => [a] -> Gen (PrimState m) -> m a
symmetricCategorical [] _ = error "symmetricCategorical: no candidates"
symmetricCategorical zs g = do
  j <- uniformR (0, length zs - 1) g
  return $ zs !! j

-- | Alias for fromIntegral.
fi :: (Integral a, Num b) => a -> b
fi = fromIntegral

max_stars_repo_name: jtobin/hnuts

max_issues_repo_name: jtobin/hnuts

max_forks_repo_name: jtobin/hnuts


text:
module Trade.Example.Correlation where

import Control.Applicative (liftA2)

import qualified Data.Set as Set

import qualified Data.Map as Map
import Data.Map (Map)

import qualified Data.Vector as Vec
import Data.Vector (Vector)

import qualified Data.ByteString.Lazy.Char8 as BSL

import Data.Time.Clock (getCurrentTime)

import qualified Statistics.Sample as Sample

import Trade.Type.BarLength (BarLength(..))

import Trade.Type.OHLC (Close(..))

import qualified Trade.Timeseries.Binance.Binance as Bin
import qualified Trade.Timeseries.Binance.Database as Bin
import qualified Trade.Timeseries.Binance.Interval as Bin
import qualified Trade.Timeseries.Binance.Symbol as Bin

import qualified Trade.Report.Heatmap as Heat

import Trade.Report.HtmlReader (render)

barLength :: BarLength
barLength = Min 1

ts :: [Bin.Symbol]
-- ts = map Bin.ETH [ Bin.ADAETH .. Bin.ZRXETH]
ts = map Bin.ETH [ Bin.BNBETH, Bin.ENJETH, Bin.HOTETH, Bin.FUELETH, Bin.XRPETH, Bin.ADAETH, Bin.EOSETH, Bin.VETETH, Bin.TRXETH, Bin.DENTETH]

getSymbols :: [Bin.Symbol] -> IO (Map Bin.Symbol (Vector Double))
getSymbols syms = do

  now <- getCurrentTime

  let req sym = Bin.RequestParams {
        Bin.baseUrl = Bin.binanceBaseUrl
        , Bin.symbol = sym
        , Bin.interval = Bin.Interval barLength
        , Bin.limit = Just 1000
        , Bin.from = Nothing
        , Bin.to = Just now
        }

      toSignal row = (Bin.toDate row, unClose (Bin.close row))

  tickers <- fmap (fmap (Vec.map toSignal)) (Bin.getTickers req syms)

  let
      tmsTickers = fmap (Set.fromList . Vec.toList . Vec.map fst) tickers
      allTimes = foldMap id tmsTickers
      tms = Map.foldr (\vs acc -> Set.intersection acc vs) allTimes tmsTickers

      res = fmap (Vec.map snd . Vec.filter ((`Set.member` tms) . fst)) tickers

 -- error (show res)
  return res

example :: IO ()
example = do
  ms <- getSymbols ts

  let xs = Map.toList ms
      f (Bin.ETH u, us) (Bin.ETH v, vs) = ((u, v), Sample.correlation (Vec.zip us vs))
      as = filter (\((x, y), _) -> x /= y) (liftA2 f xs xs)

  t <- render (Heat.heatmap 0 (fmap Just (Map.fromList as)))

  BSL.putStrLn t

max_stars_repo_name: fphh/trade

max_issues_repo_name: fphh/trade

max_forks_repo_name: fphh/trade


text:
{-# Language BangPatterns #-}
{-# Language ScopedTypeVariables #-}

module Neuralnetwork where

import           Data.Bifunctor                 ( bimap )
import           Data.Functor.Identity
import           Data.IDX
import           Data.Maybe                     ( fromMaybe )
-- import           Debug.Trace
import           Control.Monad.Random
import qualified Data.Vector.Unboxed           as V
import           Numeric.LinearAlgebra          ( Matrix
                                                , Vector
                                                , (<>)
                                                , (#>)
                                                , (><)
                                                , (|>)
                                                )
import qualified Numeric.LinearAlgebra         as LA

-- The main crux of my flaw seems to be that my weights increase without bound.
--
-- I've run into this issue _every_ single time I have written up this neural
-- network.  I've rewritten the weight updates countless times now and it's
-- starting to get to the point where my being past the deadline is a bit
-- ridiculous and if I keep it up I won't even get to the second programming
-- assignment at all.
--
-- See the delta functions and the weight update function for the most likely
-- causes of me screwing up.

-- Weights :
--   - biases = n vector of length = num of inputs to that layer
--   - nodes  = n x m where m = num of outputs to that layer
--   => m -> n layer
data Weights = W
  { wBiases :: !(Vector Double) -- length = num of inputs
  , wNodes  :: !(Matrix Double)
  } deriving Show

-- Type aliases offer no type safety but help make type signatures more
-- readable.
type Network    = [Weights]
type Input      = Vector Double
type Target     = Vector Double
type Output     = Vector Double
type Activation = Vector Double
type Prev a     = a -- Cute trick to help make type signatures more readable
type Diff a     = a

-- Here we see that we have, behold, dreaded globals.
-- However, this is fine because
--  a) I don't care
--  b) The reader pattern is overly complicated for this particular usecase
η :: Double -- Learning Rate
η = 0.1

α :: Double -- Momentum Rate
α = 0.9

-- Layer + Network generation
-- Generates a layer of random numbers between -0.05 and 0.05
randLayer :: MonadRandom m => Int -> Int -> m Weights
randLayer i o = do
  r1 :: Int <- getRandom
  r2 :: Int <- getRandom
  let bias    = (LA.randomVector r1 LA.Uniform o - 0.5) / 10
      weights = LA.uniformSample r2 o (replicate i (-0.05, 0.05))
  pure $ W bias weights

-- Generates a layer consisting of all zeroed weights
zeroLayer :: Int -> Int -> Identity Weights
zeroLayer i o = pure $ W (o |> repeat 0) (o >< i $ repeat 0)

-- A network is a list of layers; this simply steps through the list
randNet :: MonadRandom m => Int -> [Int] -> Int -> m Network
randNet i []       o = (:) <$> randLayer i o <*> pure []
randNet i (h : hs) o = (:) <$> randLayer i h <*> randNet h hs o

zeroNet :: Int -> [Int] -> Int -> Identity Network
zeroNet i []       o = (:) <$> zeroLayer i o <*> pure []
zeroNet i (h : hs) o = (:) <$> zeroLayer i h <*> zeroNet h hs o

-- Forward propagate a neural network
-- a #> b is matrix vector product where vector b is Nx1 dimensional
-- i : a 784 length vector
-- wN : (20,784) shaped matrix
runLayer :: Weights -> Input -> Activation
runLayer (W !wB !wN) !i = (wN #> i) + wB
-- runLayer (W !wB !wN) !i = trace "runLayer" $ (wN #> i) + wB

-- Gets the answer from a neural network by propagating the input forward
runNet :: Network -> Input -> Output
runNet !n !i = last $ runNet' n i
-- runNet !n !i = trace "runNet" $ last $ runNet' n i

-- Collects all intermediate activations into a list of activations
runNet' :: Network -> Input -> [Activation]
runNet' []        !i = [i]
runNet' (w : net) !i = i : runNet' net α where α = activation w i
-- runNet' (w : net) !i = trace "runNet'" $ i : runNet' net α where α = activation w i

-- Runs the neural network over a list of inputs and collects every "answer"
collectOutputs :: Network -> [Input] -> [Output]
collectOutputs net = fmap (runNet net)

logistic :: Floating a => a -> a
logistic x = 1 / (1 + exp (-x))

logistic' :: Floating a => a -> a
logistic' x = o * (1 - o) where o = logistic x

activation :: Weights -> Input -> Activation
activation !w !v = logistic $ runLayer w v

-- (Not used when actually running the neural network. The slides didn't
-- mention using a loss function at all to update the weights...)
-- Target -> Output -> Error
loss :: Target -> Output -> Double
loss = (* 0.5) ... LA.sumElements ... (^ 2) ... (-)

-- Not used anywhere in the codebase
lossO :: Floating a => a -> a -> a -> a
lossO t y α = logistic' y * (α - t)

-- If there's an error in this code, it is likely to be in the delta functions
-- or the weight updates
δo :: Target -> Output -> Vector Double
δo !t !α = α * (1 - α) * (t - α)
-- δo !t !α = trace "δo" $ α * (1 - α) * (t - α)

δh :: Weights -> Activation -> Vector Double -> Vector Double
δh (W !wB !wN) !h !δ' = h * (1 - h) * (LA.tr' wN #> δ')
-- δh (W !wB !wN) !h !δ' = trace "δh" $ h * (1 - h) * (LA.tr' wN #> δ')
-- Can't be this because it's not an inconsistent matrix error...

-- The error terms; algorithm for this and the deltas given by the slides
errors :: Network -> Target -> [Input] -> [Vector Double]
errors [!ih, !ho] !t [!αH, !αO] = [δhidden, δoutput]
-- errors [!ih, !ho] !t [!αH, !αO] = trace "errors" [δhidden, δoutput]
 where
  δoutput = δo t αO
  δhidden = δh ho αH δoutput

-- Weight deltas with momentum (only the deltas, not the actual change)
-- scalar is a function that allows me to multiply a matrix (or vector) by a single number
-- δ and i need to be the same size. δ should be 20
-- if wNΔ is (n,m) size, LA.asRow δ*i needs to be (x,m), LA.asColumn δ*i needs to be (n,x)
weightΔ :: Vector Double -> Input -> Prev (Diff Weights) -> (Weights, Diff Weights)
weightΔ !δ !i w@(W !wBΔ !wNΔ) = (w, W wΔb wΔn)
-- weightΔ !δ !i w@(W !wBΔ !wNΔ) = trace "weightΔ" (w, W wΔb wΔn)
 where
  wΔb = LA.scalar η * δ + LA.scalar α * wBΔ
  wΔn = LA.asColumn (LA.scalar η * δ * i) + LA.scalar α * wNΔ
  -- tracePrint = "δ: "      ++ show (LA.size δ)
  --           ++ " i: "     ++ show (LA.size i)
  --           ++ " wNΔ: "   ++ show (LA.size wNΔ)
  --           ++ " wBΔ: "   ++ show (LA.size wBΔ)
            -- ++ " asCol: " ++ show (LA.size $ LA.scalar η * δ * i)
            -- δ: 20 i: 784 wNΔ: (20,784) wBΔ: 20
            -- passing in wrong δ?

wgtUpd :: Weights -> Diff Weights -> Weights
wgtUpd (W wB wN) (W wB' wN') = W (wB + wB') (wN + wN')
-- wgtUpd (W wB wN) (W wB' wN') = trace "wgtUpd" $ W (wB + wB') (wN + wN')

backprop -- Given the list of prev diff weights for every layer as well
  :: (Prev Network, Network)
  -> Target
  -> Input
  -> [Prev (Diff Weights)]
  -> ([Diff Weights], Network)
backprop (net', net) !target !input w' = (wΔs, zipWith wgtUpd net wΔs)
-- backprop (net', net) !target !input w' = trace "Backprop" (wΔs, zipWith wgtUpd net wΔs)
 where
  αs  = tail $ runNet' net input -- the activations, excluding the input vector
  δs  = errors net target αs     -- length 2 lst
  wΔs = snd <$> zipWith3 weightΔ δs αs w'

-- Iterate through a list of inputs (and list of targets) and successively
-- train on those inputs.
epoch
  :: (Prev Network, Network)   -- Old network (zeros in first run), current network
  -> [Target]                  -- List of input vectors
  -> [Input]                   -- List of target vectors
  -> [Diff Weights]
  -> ([Diff Weights], Network) -- Resulting network
epoch (_    , net ) []       []       wΔs   = (wΔs, net)
epoch (!net', !net) (t : ts) (i : is) !wΔs' = epoch (net, next) ts is wΔs
-- epoch (!net', !net) (t : ts) (i : is) !wΔs' = trace "Epoch happening" $ epoch (net, next) ts is wΔs
  where (wΔs, next) = backprop (net', net) t i wΔs'

train
  :: Int -- Epochs to train for
  -> (Prev Network, Network) -- (Rest of params are identical to epoch)
  -> ([Target], [Target])
  -> ([Input], [Input])
  -> [Diff Weights]
  -> IO ()
train 0 nets (tstT, t) (tstI, i) wΔs' = do
  let final    = snd $ epoch nets t i wΔs'
      outputs  = collectOutputs final i
      oTest    = collectOutputs final tstI
      accuracy = totalError (t, outputs)
      tstAcc   = totalError (tstT, oTest)
  print $ "Epoch 50 accuracy: " ++ show accuracy
  print $ "Epoch 50 test accuracy: " ++ show tstAcc
  -- print out stuff
train n (net', net) t'@(tstT, t) i'@(tstI, i) wΔs' = do
  let (wΔs, next) = epoch (net', net) t i wΔs'
      outputs     = collectOutputs next i
      oTest       = collectOutputs next tstI
      accuracy    = totalError (t, outputs)
      tstAcc      = totalError (tstT, oTest)
  print $ "Epoch " ++ show (50-n) ++ " accuracy: "      ++ show accuracy
  print $ "Epoch " ++ show (50-n) ++ " test accuracy: " ++ show tstAcc
  train (n - 1) (net, next) t' i' wΔs

-- Simple debugging function that tells me the shape of a particular network.
-- Useful for ensuring my matrices are the right shape
sizeOf :: Network -> [(Int, (Int, Int))]
sizeOf [W hB hN, W oB oN] = [(LA.size hB, LA.size hN), (LA.size oB, LA.size oN)]

totalError :: ([Target], [Output]) -> Double
totalError (!ts, !os) = sum (zipWith loss ts os) / fromIntegral (length ts)

main :: IO ()
main = do
  -- Yes, this is slightly ugly...
  !file <- fromMaybe (error "file decoding failed")
    <$> decodeIDXFile "/home/jaredweakly/Documents/Classes/CS445/train-images-idx3-ubyte"
  !labels <- fromMaybe (error "label decoding failed") <$> decodeIDXLabelsFile
    "/home/jaredweakly/Documents/Classes/CS445/train-labels-idx1-ubyte"
  !testFile <- fromMaybe (error "test file decoding failed")
    <$> decodeIDXFile "/home/jaredweakly/Documents/Classes/CS445/t10k-images-idx3-ubyte"
  !testLabels <- fromMaybe (error "test label decoding failed") <$> decodeIDXLabelsFile
    "/home/jaredweakly/Documents/Classes/CS445/t10k-labels-idx1-ubyte"

  -- mnist : 60_000 [(target, Vector 784 Double)]
  let !mnist = fromMaybe (error "labeling failed") $ labeledDoubleData labels file
      !mnist' =
        fromMaybe (error "labeling failed") $ labeledDoubleData testLabels testFile
      !initN20   = runIdentity $ zeroNet 784 [20] 10
      !initN50   = runIdentity $ zeroNet 784 [50] 10
      !initN100  = runIdentity $ zeroNet 784 [100] 10
      !testset = mnist' ++ mnist

  print "Loaded up data"
  -- inputs : (60000><785) Matrix Double
  -- targets : Vector 10 Double. 0.1 everywhere but target which is 0.9
  let mkSparseVectors [] = [] :: [Vector Double]
      mkSparseVectors (x : xs) =
        LA.assoc 10 0.1 [(x, 0.9)] : mkSparseVectors xs :: [Vector Double]

      normalize :: V.Vector Double -> Vector Double
      normalize            = (/ 255) . V.convert

      !(!targets, !inputs) = bimap mkSparseVectors id . unzip $ normalize <$$> mnist
      !(!tstT   , !tstI  ) = bimap mkSparseVectors id . unzip $ normalize <$$> testset

  !weights20  <- randNet 784 [20] 10
  !weights50  <- randNet 784 [50] 10
  !weights100 <- randNet 784 [100] 10

  print "Beginning training with n=20"
  train 50 (initN20, weights20) (tstT, targets) (tstI, inputs) initN20

  print "Beginning training with n=50"
  train 50 (initN50, weights50) (tstT, targets) (tstI, inputs) initN50

  print "Beginning training with n=100"
  train 100 (initN100, weights100) (tstT, targets) (tstI, inputs) initN100

-- Helper functions I was too lazy to import
-- These are not essential to understanding any of the actual code.

-- Lets me write pointfree code using functions that takes two arguments
-- instead of just one
infixr 9 ...
(...) :: (b -> c) -> (a1 -> a2 -> b) -> a1 -> a2 -> c
(...) = (.) . (.)

-- Lets me map over something inside another structure.
-- eg: (+1) <$$> [[a]] would map over all of the inner lists
infixl 4 <$$>
(<$$>) :: (Functor f2, Functor f1) => (a -> b) -> f1 (f2 a) -> f1 (f2 b)
(<$$>) = fmap fmap fmap

both :: (a -> b) -> (a, a) -> (b, b)
both f (a, a') = (f a, f a')

max_stars_repo_name: jared-w/Haskell-MNIST

max_issues_repo_name: jared-w/Haskell-MNIST

max_forks_repo_name: jared-w/Haskell-MNIST


text:
module Main where

import Data.Complex

main :: IO ()
main = mapM_ (putStrLn . line) ys
  where
    line y = map (\z -> if z < escape && not (isNaN z) then ' ' else '#') [magnitude $ man (x :+ y) | x <- xs]
    man = mandelbrot iters (seed :+ 0)
    iters = 50
    seed = 0.0099
    escape = 2
    bresx = 50
    bresy = 30
    scale = 1
    xpos = -0.5
    ypos = 0
    xs = map (\x -> x*scale + xpos) [x/bresx | x <- [-bresx..bresx]]
    ys = map (\y -> y*scale + ypos) [y/bresy | y <- [-bresy..bresy]]

-- | n Mandelbrot iterations starting with initial value z and constant c.
mandelbrot :: Int -> Complex Double -> Complex Double -> Complex Double
mandelbrot n z c
  | n <= 0 = z
  | otherwise = mandelbrot (pred n) (z * z + c) c

max_stars_repo_name: Kludgy/mandeltoy-hs

max_issues_repo_name: Kludgy/mandeltoy-hs

max_forks_repo_name: Kludgy/mandeltoy-hs


text:
{-# LANGUAGE DeriveDataTypeable, PatternGuards #-}
module Physics.Quaternion
  ( Quaternion(..)
  , Complicated(..)
  , Hamiltonian(..)
  , slerp
  , asinq
  , acosq
  , atanq
  , asinhq
  , acoshq
  , atanhq
  , absi
  , pow
  , rotate
  ) where

import Control.Applicative
import Control.Lens
import Control.Lens.Rep
import Data.Complex (Complex((:+)))
import Data.Data
import Data.Distributive
import Data.Foldable
import Data.Monoid
import Data.Traversable
import Physics.Epsilon
import Physics.Involutive
import Physics.Metric
import Physics.V3
import Physics.Vector
import Prelude hiding (any)

data Quaternion a = Quaternion a a a a deriving (Eq,Ord,Read,Show,Data,Typeable)

instance Functor Quaternion where
  fmap f (Quaternion e i j k) = Quaternion (f e) (f i) (f j) (f k)
  a <$ _ = Quaternion a a a a

instance Applicative Quaternion where
  pure a = Quaternion a a a a
  Quaternion f g h i <*> Quaternion a b c d = Quaternion (f a) (g b) (h c) (i d)

instance Monad Quaternion where
  return = pure
  (>>=) = bindRep

instance Representable Quaternion where
  rep f = Quaternion (f e) (f i) (f j) (f k)

instance Foldable Quaternion where
  foldMap f (Quaternion e i j k) = f e `mappend` f i `mappend` f j `mappend` f k
  foldr f z (Quaternion e i j k) = f e (f i (f j (f k z)))

instance Traversable Quaternion where
  traverse f (Quaternion e i j k) = Quaternion <$> f e <*> f i <*> f j <*> f k

instance RealFloat a => Num (Quaternion a) where
  {-# SPECIALIZE instance Num (Quaternion Float) #-}
  {-# SPECIALIZE instance Num (Quaternion Double) #-}
  (+) = liftA2 (+)
  (-) = liftA2 (-)
  negate = fmap negate
  Quaternion a1 b1 c1 d1 * Quaternion a2 b2 c2 d2 = Quaternion
    (a1*a2 - b1*b2 - c1*c2 - d1*d2)
    (a1*b2 + b1*a2 + c1*d2 - d1*c2)
    (a1*c2 - b1*d2 + c1*a2 + d1*b2)
    (a1*d2 + b1*c2 - c1*b2 + d1*a2)
  fromInteger x = Quaternion (fromInteger x) 0 0 0
  abs z = Quaternion (norm z) 0 0 0
  signum q@(Quaternion e i j k)
    | m == 0.0 = q
    | not (isInfinite m || isNaN m) = q ^/ sqrt m
    | any isNaN q = qNaN
    | not (ii || ij || ik) = Quaternion 1 0 0 0
    | not (ie || ij || ik) = Quaternion 0 1 0 0
    | not (ie || ii || ik) = Quaternion 0 0 1 0
    | not (ie || ii || ij) = Quaternion 0 0 0 1
    | otherwise = qNaN
    where
      m = quadrance q
      ie = isInfinite e
      ii = isInfinite i
      ij = isInfinite j
      ik = isInfinite k

  -- abs    = error "Quaternion.abs: use norm"
  -- signum = error "Quaternion.signum: use signorm"

qNaN :: RealFloat a => Quaternion a
qNaN = Quaternion fNaN fNaN fNaN fNaN where fNaN = 0/0

-- {-# RULES "abs/norm" abs x = Quaternion (norm x) 0 0 0 #-}
-- {-# RULES "signum/signorm" signum = signorm #-}

-- this will attempt to rewrite calls to abs to use norm intead when it is available.

instance RealFloat a => Fractional (Quaternion a) where
  {-# SPECIALIZE instance Fractional (Quaternion Float) #-}
  {-# SPECIALIZE instance Fractional (Quaternion Double) #-}
  Quaternion q0 q1 q2 q3 / Quaternion r0 r1 r2 r3 = Quaternion (r0*q0+r1*q1+r2*q2+r3*q3)
                                                               (r0*q1-r1*q0-r2*q3+r3*q2)
                                                               (r0*q2+r1*q3-r2*q0-r3*q1)
                                                               (r0*q3-r1*q2+r2*q1-r3*q0)
                                                 ^/ (r0*r0 + r1*r1 + r2*r2 + r3*r3)
  recip q = q ^/ quadrance q
  fromRational x = Quaternion (fromRational x) 0 0 0

instance Metric Quaternion where
  Quaternion e i j k `dot` Quaternion e' i' j' k' = e*e' + i*i' + j*j' + k*k'

class Complicated t where
  e :: Functor f => (a -> f a) -> t a -> f (t a)
  i :: Functor f => (a -> f a) -> t a -> f (t a)

instance Complicated Complex where
  e f (a :+ b) = (:+ b) <$> f a
  i f (a :+ b) = (a :+) <$> f b

instance Complicated Quaternion where
  e f (Quaternion a b c d) = (\a' -> Quaternion a' b c d) <$> f a
  i f (Quaternion a b c d) = (\b' -> Quaternion a b' c d) <$> f b

class Complicated t => Hamiltonian t where
  j :: Functor f => (a -> f a) -> t a -> f (t a)
  k :: Functor f => (a -> f a) -> t a -> f (t a)
  ijk :: Functor f => (V3 a -> f (V3 a)) -> t a -> f (t a)

instance Hamiltonian Quaternion where
  j f (Quaternion a b c d) = (\c' -> Quaternion a b c' d) <$> f c
  k f (Quaternion a b c d) = Quaternion a b c <$> f d

  ijk f (Quaternion a b c d) = (\(V3 b' c' d') -> Quaternion a b' c' d') <$> f (V3 b c d)

instance Distributive Quaternion where
  distribute f = Quaternion (fmap (^.e) f) (fmap (^.i) f) (fmap (^.j) f) (fmap (^.k) f)

instance (Involutive a, Num a) => Involutive (Quaternion a) where
  conjugate (Quaternion e i j k) = Quaternion (conjugate e) (-i) (-j) (-k)

reimagine :: RealFloat a => a -> a -> Quaternion a -> Quaternion a
reimagine r s (Quaternion _ i j k)
  | isNaN s || isInfinite s = Quaternion r
    (if i /= 0 then i * s else 0)
    (if j /= 0 then j * s else 0)
    (if k /= 0 then k * s else 0)
  | otherwise = Quaternion r (i * s)  (j * s) (k * s)

-- | quadrance of the imaginary component
qi :: Num a => Quaternion a -> a
qi (Quaternion _ i j k) = i*i + j*j + k*k

absi :: Floating a => Quaternion a -> a
absi = sqrt . qi

pow :: RealFloat a => Quaternion a -> a -> Quaternion a
pow q t = exp (t *^ log q)

-- ehh..
instance RealFloat a => Floating (Quaternion a) where
  {-# SPECIALIZE instance Floating (Quaternion Float) #-}
  {-# SPECIALIZE instance Floating (Quaternion Double) #-}
  pi = Quaternion pi 0 0 0
  exp q@(Quaternion e i j k)
    | qiq == 0 = Quaternion (exp e) i j k
    | ai <- sqrt qiq, ee <- exp e = reimagine (ee * cos ai) (ee * (sin ai / ai)) q
    where qiq = qi q
  log q@(Quaternion e i j k)
    | qiq == 0 = if e >= 0 then Quaternion (log e) i j k else Quaternion (log (negate e)) pi j k -- mmm, pi
    | ai <- sqrt qiq, m <- sqrt (e*e + qiq) = reimagine (log m) (atan2 m e / ai) q
    where qiq = qi q
  x ** y = exp (y * log x)
  sqrt q@(Quaternion e i j k)
    | m   == 0 = q
    | qiq == 0 = if e > 0 then Quaternion (sqrt e) 0 0 0 else Quaternion 0 (sqrt (negate e)) 0 0
    | im <- sqrt (0.5*(m-e)) / sqrt qiq = Quaternion (0.5*(m+e)) (i*im) (j*im) (k*im)
    where qiq = qi q
          m = sqrt (e*e + qiq)
  cos q@(Quaternion e i j k)
    | qiq == 0 = Quaternion (cos e) i j k
    | ai <- sqrt qiq = reimagine (cos e * cosh ai) (- sin e * (sinh ai / ai)) q
    where qiq = qi q
  sin q@(Quaternion e i j k)
    | qiq == 0 = Quaternion (sin e) i j k
    | ai <- sqrt qiq = reimagine (sin e * cosh ai) (cos e * (sinh ai / ai)) q
    where qiq = qi q
  tan q@(Quaternion e i j k)
    | qiq == 0 = Quaternion (tan e) i j k
    | ai <- sqrt qiq, ce <- cos e, sai <- sinh ai, d <- ce*ce + sai*sai =
      reimagine (ce * sin e / d) (cosh ai * (sai / ai) / d) q
    where qiq = qi q
  sinh q@(Quaternion e i j k)
    | qiq == 0 = Quaternion (sinh e) i j k
    | ai <- sqrt qiq = reimagine (sinh e * cos ai) (cosh e * (sin ai / ai)) q
    where qiq = qi q
  cosh q@(Quaternion e i j k)
    | qiq == 0 = Quaternion (cosh e) i j k
    | ai <- sqrt qiq = reimagine (cosh e * cos ai) ((sinh e * sin ai) / ai) q
    where qiq = qi q
  tanh q@(Quaternion e i j k)
    | qiq == 0 = Quaternion (tanh e) i j k
    | ai <- sqrt qiq, se <- sinh e, cai <- cos ai, d <- se*se + cai*cai =
      reimagine ((cosh e * se) / d) ((cai * (sin ai / ai)) / d) q
    where qiq = qi q

  asin q = cut asin q
  acos q = cut acos q
  atan q = cut atan q

  asinh q = cut asinh q
  acosh q = cut acosh q
  atanh q = cut atanh q

cut :: RealFloat a => (Complex a -> Complex a) -> Quaternion a -> Quaternion a
cut f q@(Quaternion e _ j k)
  | qiq == 0 = Quaternion a b j k
  | otherwise = reimagine a (b / ai) q
  where qiq = qi q
        ai = sqrt qiq
        a :+ b = f (e :+ ai)

cutWith :: RealFloat a => Complex a -> Quaternion a -> Quaternion a
cutWith (r :+ im) q@(Quaternion e i j k)
  | e /= 0 || qiq == 0 || isNaN qiq || isInfinite qiq = error "bad cut"
  | s <- im / sqrt qiq = Quaternion r (i*s) (j*s) (k*s)
  where qiq = qi q

asinq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
asinq q@(Quaternion e _ _ _) u
  | qiq /= 0.0 || e >= -1 && e <= 1 = asin q
  | otherwise = cutWith (asin (e :+ sqrt qiq)) u
  where qiq = qi q

acosq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
acosq q@(Quaternion e _ _ _) u
  | qiq /= 0.0 || e >= -1 && e <= 1 = acos q
  | otherwise = cutWith (acos (e :+ sqrt qiq)) u
  where qiq = qi q

atanq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
atanq q@(Quaternion e _ _ _) u
  | e /= 0.0 || qiq >= -1 && qiq <= 1 = atan q
  | otherwise = cutWith (atan (e :+ sqrt qiq)) u
  where qiq = qi q

asinhq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
asinhq q@(Quaternion e _ _ _) u
  | e /= 0.0 || qiq >= -1 && qiq <= 1 = asinh q
  | otherwise = cutWith (asinh (e :+ sqrt qiq)) u
  where qiq = qi q

acoshq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
acoshq q@(Quaternion e _ _ _) u
  | qiq /= 0.0 || e >= 1 = asinh q
  | otherwise = cutWith (acosh (e :+ sqrt qiq)) u
  where qiq = qi q

atanhq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
atanhq q@(Quaternion e _ _ _) u
  | qiq /= 0.0 || e > -1 && e < 1 = atanh q
  | otherwise = cutWith (atanh (e :+ sqrt qiq)) u
  where qiq = qi q

slerp :: RealFloat a => Quaternion a -> Quaternion a -> a -> Quaternion a
slerp q p t
  | 1.0 - cosphi < 1e-8 = q
  | phi <- acos cosphi, r <- recip (sin phi)
  = (sin ((1-t)*phi)*r *^ q ^+^ f (sin (t*phi)*r) *^ p) ^/ sin phi
  where
   dqp = dot q p
   (cosphi, f) = if dqp < 0 then (-dqp, negate) else (dqp, id)
{-# SPECIALIZE slerp :: Quaternion Float -> Quaternion Float -> Float -> Quaternion Float #-}
{-# SPECIALIZE slerp :: Quaternion Double -> Quaternion Double -> Double -> Quaternion Double #-}

--slerp :: RealFloat a => Quaternion a -> Quaternion a -> a -> Quaternion a
--slerp q0 q1 = let q10 = q1 / q0 in \t -> pow q10 t * q0

rotate :: Num a => Quaternion a -> V3 a -> V3 a
rotate (Quaternion a' b c d) (V3 x y z) = V3
  (2*((t8+t10)*x+(t6- t4)*y+(t3+t7)*z)+x)
  (2*((t4+ t6)*y+(t5+t10)*y+(t9-t2)*z)+y)
  (2*((t7- t3)*z+(t2+ t9)*z+(t5+t8)*z)+z)
  where
    a = -a'
    t2 = a*b
    t3 = a*c
    t4 = a*d
    t5 = -b*b
    t6 = b*c
    t7 = b*d
    t8 = -c*c
    t9 = c*d
    t10 = -d*d
{-# SPECIALIZE rotate :: Quaternion Float -> V3 Float -> V3 Float #-}
{-# SPECIALIZE rotate :: Quaternion Double -> V3 Double -> V3 Double #-}

instance (RealFloat a, Epsilon a) => Epsilon (Quaternion a) where
  nearZero = nearZero . quadrance

max_stars_repo_name: ekmett/physics

max_issues_repo_name: ekmett/physics

max_forks_repo_name: ekmett/physics