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

module FEECa.Internal.Form (

  -- * Generic form types
    Dim, Form (Form, arity, dimVec, terms), split

  -- * Predefined primitive constructors
  , zeroForm, nullForm, oneForm

  -- * Form operations
  , apply, refine, refineBasis, inner, contract, trace

  ) where

import            Data.List ( intersect, elemIndex )
import            Data.Maybe ( fromJust )
import qualified  Numeric.LinearAlgebra.HMatrix as M
-- import qualified Numeric.LinearAlgebra.Data as M

import            FEECa.Utility.Combinatorics
import            FEECa.Utility.Discrete
import qualified  FEECa.Utility.Print           as P  ( Pretty(..), printForm,
                                                        text, (<>), (<+>), int )
import            FEECa.Utility.Utility               ( pairM, sumV, expSign, sign )

import            FEECa.Internal.Spaces     hiding    ( inner )
import qualified  FEECa.Internal.Spaces         as S  ( inner )

-- * General form: does not depend on the underlying vector space it works on
--   in any way.

type Dim = Int
type Idx = Int
type Prd = [Idx] -- a product of projections / indexed 1-forms
-- type LinComb f = [(f, Prd)]

-- | Bilinear, alternating forms over vectorspaces
data Form f =  -- we lose dependency on the type of vector!
    Form  { arity   :: Dim            -- ^ For complete evaluation
          , dimVec  :: Dim            -- ^ Of the underlying vector space
          , terms   :: [(f, Prd)] }   -- ^ List of terms of (coeff,wedge)'s
  deriving Eq

split :: (Module w, Scalar w ~ v) => Form w -> ([w], [Form v])
split (Form k n cs) = unzip $ map split' cs
  where split' (a,b) = (a, Form k n [(mulId, b)])

-- terms [(17, [1,2]), (38, [1,3])] = 17*dx1/\dx2 + 38*dx1/\dx3

-- | Invariant for terms of a defined form: all terms have the same arity
--   ie: each is the result of the exterior product of the same number of 1-forms
termsInv :: [(f, Prd)] -> Bool
termsInv []          = True
termsInv ((_,xs):ys) = all (\(_,xs') -> length xs == length xs') ys

-- NB: because of ((,) t) 's functorial nature, maybe it could make sense to
--   rearrange our terms so as to have them be (inds,coeff) like they used to be?
-- XXX: change so as to inspect result (in case f zeroes out a coeff?)
instance Functor Form where
  fmap f (Form k n cs) = Form k n (map (pairM f id) cs)

instance Show f => Show (Form f) where
  show (Form k n cs) = show k ++ "-form in " ++ show n ++ " dimensions: " ++
                        show cs

instance P.Pretty f => P.Pretty (Form f) where
  pPrint (Form k n cs) =
    P.int k P.<> P.text "-form in" P.<+> P.int n P.<+> P.text "dimensions:"
    P.<+> P.printForm "dx" "0" P.pPrint cs

-- NB: will be (i -> i -> Ordering) once we normalise all products to be in
--      their normal form - an increasing list
combineWithBy :: (a -> a -> a) -> (i -> i -> Bool)
              -> [(a, i)] -> [(a, i)] -> [(a,i)]
combineWithBy f p = foldl ins
  where ins [] x = [x]
        ins (y:ys) x | p (snd x) (snd y) = (fst x `f` fst y, snd y) : ys
                     | otherwise         = y : ins ys x

-- XXX: have combinator
--        aggregate f p cs1 cs2
--      where f :: a -> a -> a (or b -> c)
--            p :: MultiIndex -> MultiIndex -> Bool
--            cs1, cs2 :: [(a, MultiIndex)]
--      such that combines elements from cs2 with those in cs1 (sequentially,
--      no repeat, or filtering those appropriate) with p defining the "quotient"
--      class, f defining the transformation on coeff's
--      OPT: also a MultiIndex transformation? or a separte combinator for it?

-- | Sum of forms
-- Shall we do: permutation simplification/identification
(+++) :: Ring f => Form f -> Form f -> Form f
omega +++ eta
    | degNEq omega eta = errForm "(+++)" BiDegEq
    | spaNEq omega eta = errForm "(+++)" BiSpaEq
    | otherwise = Form (arity eta) (dimVec eta)
                       (step (terms omega) (terms eta))
  where step [] ys = ys
        step xs [] = xs
        step (x:xs) (y:ys)
          | snd x == snd y = let z = add (fst x) (fst y) in
              if fst x /= (addInv . fst) y then (z, snd x) : step xs ys
                            else step xs ys
          | snd x < snd y  = x : step xs (y:ys)
          | otherwise      = y : step (x:xs) ys

-- | Scaling of forms
(***) :: Ring f => f -> Form f -> Form f
(***) a | a == addId = \(Form k n _) -> zeroForm k n
        | otherwise  = fmap (mul a)

-- | (Exterior) Product of forms
(//\\) :: Ring f => Form f -> Form f -> Form f
omega //\\ eta
    | spaNEq omega eta = errForm "(//\\\\)" BiSpaEq
    | otherwise = Form (arity omega + arity eta) (dimVec eta)
                       (concatMap (\d -> map (`combine` d) (dxs d)) (terms eta))
  where dxs     (_,ys) = filter (null . intersect ys . snd) (terms omega)
        combine (a,xs) = pairM (mul a) (xs++)

-- | Forms over a 'Ring' form a 'Module'.
instance Ring f => Module (Form f) where
  type Scalar (Form f) = f
  addV = (+++)
  sclV = (***)

-- | Forms over a 'Field' form a 'VectorSpace'.
instance Field f => VectorSpace (Form f)

-- | For 'Form's defined over a 'Ring' we associate an 'Algebra': the exterior
-- algebra.
-- This instance is valid this way since we do not have a restrictive typing
-- for forms and hence addition is blind to arity *type-wise* - however,
-- runtime errors will take place if invalid addition is attempted.
instance Ring f => Algebra (Form f) where
  addA = addV
  (/\) = (//\\)
  sclA = sclV

-- | Basic abstract 1-form
oneForm :: Ring f => Dim -> Dim -> Form f
oneForm i n | i < 0 || i > n  = errForm "oneForm" MoProjBd
            | otherwise       = Form 1 n [ (mulId,[i]) ]

-- TODO: shall we have something special for these? no need to state dimension
-- n since they will be constantly zero anyway

-- | The (normalised) == 0 form
zeroForm :: Dim -> Dim -> Form f
zeroForm k n = Form k n []

-- | The k-arity == 0 form
nullForm :: Dim -> f -> Form f
nullForm n f = Form 0 n [(f, [])]

-- Necesitamos una función de pinchado
--  y así pinchar las consecutivas componentes
-- If the function was actually a field, this part would be simplified
contract :: (Ring f, Module v, Dimensioned v)
         => (Idx -> v -> f) -> Form f -> v -> Form f
contract proj omega v
    | vecNEq omega v = errForm "contract" MoVecEq
    | otherwise      = foldl (+++) (zeroForm k n) $
        map (\c -> Form k n $ map (pinchado c) [1..arity omega])
            (terms omega)
        {- concatMap (\c -> map (pinchado c) [1..arity omega]) -}
  where k = max 0 (arity omega - 1)
        n = dimVec omega
        pinchado (f,[]) _ = (f, []) -- error ??
        pinchado (f,ds) i = let (ds1,j:ds2) = splitAt (i-1) ds in
                              (expSign i $ mul f (proj j v), ds1 ++ ds2)
  {- TODO:  optimise
            error handling: proj indexing beyond dimension of v -}

-- list_to_index :: Integral a => a -> [a] -> a
-- list_to_index n (l:ls)

-- list_to_index' :: Integral a => a -> a -> [a] -> a
-- list_to_index' acc n (l:ls) = list_to_index' ((acc * n) + l) n ls
-- list_to_index' acc   _ []     = acc

-- make_lookup  :: (Ring r, EuclideanSpace v, Scalar v ~ r)
--              => Int -> Int -> [v] -> [[v]] -> Array [r]
-- make_lookup n k ds vvs
apply :: (EuclideanSpace v, Ring w, Module w, Scalar v ~ Scalar w)
      => [v] -> [v] -> Form w -> w
apply ds vs (Form k _ cs) = foldl addV addId (map (apply' k ds vs) cs)

apply' :: (EuclideanSpace v, Module w, Scalar v ~ Scalar w)
       => Int -> [v] -> [v] -> (w,[Int]) -> w
apply' _ _  _  (p, []) = p
apply' k ds vs (p, cs) = sclV c p
  where projections    = [toDouble $ dot (ds !! i) v | v <- vs, i <- cs]
        c              = fromDouble $ M.det $ M.matrix k projections

refineBasis :: (EuclideanSpace v, Scalar v ~ r)
            => [v] -> [[v]] -> [[r]]
refineBasis ds vvs = map (map (fromDouble . M.det)) submatrices
  where projections = [[[toDouble $ dot d v | v <- vs] | d <- ds] | vs <- vvs]
        matrices    = map (kSublists k) projections
        submatrices = map (map (M.matrix k . concat)) matrices
        k           = length (head vvs)

-- | Run function for 'Form's: given (an appropriate number of) vector arguments
--   and a 1-form basis (given as a basis-element indexing function 'proj'), it
--   evaluates the form on those arguments
refine :: (Ring w, Module w, Module v, Scalar v ~ Scalar w)
       => (Idx -> v -> Scalar w)  -- ^ The definition for the projection function
                                  --   for the specific vector space
       -> Form w
       -> [v] -> w
refine proj (Form _ _ cs) vs = {-#SCC "Form.refine" #-} sumV (map (($ vs) . formify proj) cs')
  where cs' | null cs   = [(addId,[])]
            | otherwise = cs
-- TODO: capture inconsistency between k and length vs here??
-- ALSO: 0-forms... not evaluating correctly now! Cfr: formify does not accept
--    empty cs
-- XXX: for now proj should take care of the error... change later when settled

-- | Helper function in evaluation: given a 1-form basis, converts a single
--   'Form' term into an actual function on vectors
formify :: (Ring w, Module w, Module v, Scalar w ~ Scalar v)
        => (i -> v -> Scalar v) -> (w,[i]) -> [v] -> w
formify _    (s, [])   _  = s
formify proj (s, i:is) vs
    | null is   = {-#SCC "Form.formify"  #-} sclV (proj i (head vs)) s
    | otherwise = {-#SCC "Form.formifyR" #-}
        foldl addV addId
              (map (\(w,e) -> sclV
                                (mul (sign (w,e)) ((proj i . head) (pick' w vs)))
                                (formify proj (s,is) (pick' e vs)))
                   (permutationPairs (length is + 1) 1 (length is)))
  where pick' ns = pick (differences ns)

-- We need a basis here
inner :: (InnerProductSpace w, EuclideanSpace v, Scalar w ~ Scalar v)
      => (Idx -> v -> Scalar w)  -- ^ Projection function in the specific vector space
      -> Form w -> Form w -> Scalar w
inner proj omega eta
    | degNEq omega eta = errForm "inner" BiDegEq -- TODO (??)
    | otherwise = foldl
          (flip $ \vs -> add (S.inner (app omega vs) (app eta vs)))
          addId
          (map pick' (permutations n (arity omega)))
  where pick' is  = pick (differences is) (map (unitVector n) [0..n-1])
        app       = refine proj
        n         = dimVec omega

trace :: [Int] -> Form w -> Form w
trace sigma (Form k n ts)
    | k' < k    = zeroForm k' n
    | otherwise = Form k n (map (pairM id restrict') ts')
  where ts'       = filter (is_in_range' sigma . snd) ts
        restrict' = map (fromJust . (`elemIndex` sigma))
        k'        = length sigma
        -- is_in_range' :: [Int] -> [Int] -> Bool
        is_in_range' sigma = all (`elem` sigma)

-- | Checks arity equality
degNEq :: Form f -> Form f -> Bool
degNEq omega eta = arity omega /= arity eta

-- | Checks combined arity bound
degNBd :: Form f -> Form f -> Bool
degNBd  omega eta = (arity omega + arity eta) <= dimVec omega

-- | Checks compatible underlying vector space dimensions between forms
spaNEq :: Form f -> Form f -> Bool
spaNEq omega eta = dimVec omega /= dimVec eta

-- | Checks compatible underlying vector space dimensions between a form and a
-- 'Dimensioned' type value
vecNEq :: Dimensioned v => Form f -> v -> Bool
vecNEq omega v = dimVec omega /= dim v

errForm :: String -> FormMust -> t
errForm callee obligation = error $ "Form." ++ callee ++
                                    ": forms must " ++ show obligation

-- | Kinds of enforcements to the definitions and operations between/for 'Form'
data FormMust = BiDegEq | BiDegBd | BiSpaEq | MoProjBd | MoVecEq

instance Show FormMust where
  show BiDegEq  = "be of the same degree"
  show BiDegBd  = "have joint degree bounded by the working vector space dimension"
  show BiSpaEq  = "act on the same vector space"
  show MoProjBd = "project components of the underlying vector space"
  show MoVecEq  = "act on vectors of the working vectors space"

max_stars_repo_name: Airini/FEECa

max_issues_repo_name: Airini/FEECa

max_forks_repo_name: Airini/FEECa


text:
module AlgorithmsSpec where

import qualified Statistics.Matrix as M
import qualified Statistics.Matrix.Types as T
import qualified Statistics.Matrix.Algorithms as Alg
import qualified Statistics.Matrix.Function as Func
import qualified Statistics.Matrix.Mutable as Mut

import Utils

qrDecompositionInvariant :: T.Matrix -> Bool
qrDecompositionInvariant mat =
  (M.multiply (fst res) (snd res)) == mat
  where
    res = Alg.qr mat

qrFirstOrthoInvariant :: T.Matrix -> Bool
qrFirstOrthoInvariant mat = Utils.isOrtho (fst res)
  where
    res = Alg.qr mat

qrSecondTriInvariant :: T.Matrix -> Bool
qrSecondTriInvariant mat = Utils.isUpperTri (snd res)
  where
    res = Alg.qr mat

ext: hs

max_stars_repo_name: Magalame/dh-core

max_issues_repo_name: Magalame/dh-core

max_forks_repo_name: Magalame/dh-core


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

import Control.Exception
import Data.Complex
import Data.Word
import Foreign.C.Types
import GHC.Int
import Test.Hspec

import ArrayFire

spec :: Spec
spec =
  describe "Image tests" $ do
    it "Should test if Image I/O is available" $ do
      isImageIOAvailable `shouldReturn` True

max_stars_repo_name: danielkroeni/arrayfire-haskell

max_issues_repo_name: danielkroeni/arrayfire-haskell

max_forks_repo_name: danielkroeni/arrayfire-haskell


text:
module Noether.Algebra.Actions.Compatible where

import           Data.Complex

import           Noether.Lemmata.Prelude
import           Noether.Lemmata.TypeFu

import           Noether.Algebra.Single
import           Noether.Algebra.Tags

import           Noether.Algebra.Actions.Acts

{-| A strategy-parameterized typeclass for a compatible action, where compatibility
    is defined in the group action sense.

    A compatible action satisfies
    a `act` (a' `act` b) = (a `op` a') `act` b
-}

class CompatibleK (lr :: Side) (op :: k1) (act :: k2) a b (s :: CompatibleE)

data CompatibleE = Compatible_Acts_Semigroup
  { compatible_actor           :: Type
  , compatible_action          :: ActsE
  , compatible_actor_semigroup :: SemigroupE
  }

type family CompatibleS (lr :: Side) (op :: k1) (act :: k2) (a :: Type) (b :: Type) = (r :: CompatibleE)

instance (ActsK lr act a b za, SemigroupK op a zs) =>
         CompatibleK lr op act a b (Compatible_Acts_Semigroup a za zs)

max_stars_repo_name: evertedsphere/noether

max_issues_repo_name: evertedsphere/noether

max_forks_repo_name: evertedsphere/noether


text:
{-# LANGUAGE CPP #-}

module Year2021.Day15 where

import Numeric.LinearAlgebra
import Util
import Safe
import Linear.V2
import System.FilePath
import Data.Bifunctor

part1, part2 :: Int
part1 = 1
part2 = 5

main = readFile (replaceExtension __FILE__ ".in") >>= \input ->
  let grid = fromLists . map (map (read . pure)) $ lines input :: Matrix Z
      grid' = cmod 9 (fromBlocks [[grid + scalar (r + c) | c <- [0..4]] | r <- [0..4]] - 1) + 1
      neigh (a, b) = [ (grid' ! x ! y, (x, y))
        | p@(V2 x y) <- (+ V2 a b) <$> adjacent
        , inBounds 0 (uncurry V2 (size grid') - 1) p ]
   in print . fst . findJust ((== size grid') . bimap succ succ . head . snd)
        $ dijkstra neigh (0, 0)


text:
{-# LANGUAGE GADTs, StandaloneDeriving, KindSignatures, FlexibleContexts, UndecidableInstances, FlexibleInstances #-}
module Data.Mat.Data where

import Control.Monad
import Control.Monad.Trans
import Data.Complex
import Data.Functor
import Data.List
import Foreign
import Foreign.C

import Data.Mat.MemParser
import Data.Mat.Raw

-- | Representation of a MAT data type.
--
-- The first argument determines whether this is a type pattern or an
-- actual type.
--
-- The second argument is the "underlying" datatype.
--
-- The third argument is the Haskell result type.
--
-- TODO: Extend the documentation.
data MatType :: (* -> *) -> * -> * -> * where
  -- Base types
  MatDouble  :: MatType d Double Double  -- TODO: should this be CDouble?
  MatInt64   :: MatType d Int64  Int64   -- TODO: should this be CInt64?
  MatBool    :: MatType d Word8  Bool    -- TODO: should this be CUInt8?
  MatChar    :: MatType d CUChar Char

  -- Dynamic type
  MatCell    :: MatType d (Ptr CMatVar) MatData
    -- ^ Represents a single cell, i.e., a single value of arbitrary type.

  -- Compound types
  MatComplex :: MatType d a b -> MatType d (CMatComplexSplit a) (Complex b)
  MatStruct  :: d [String] -> MatType d (Ptr CMatVar) MatStruct
    -- ^ Represents a struct array with the given number of fields.
    -- TODO: The representation should contain the names of the fields.

  MatArray   :: MatType d a b -> d Int -> MatType d a [b]
    -- ^ Represents an array with the given element type and number of
    -- elements.

  -- Other
  MatUnknown :: MatioClass -> MatioType -> MatType d (Ptr CData) ()
    -- ^ Represents a value we do not know how to interpret properly as
    -- a Haskell type. We store the class and type of the value.

showMatType :: MatType Known a b -> String
showMatType MatDouble              = "Double"
showMatType MatInt64               = "Int64"
showMatType MatBool                = "Bool"
showMatType MatChar                = "Char"
showMatType MatCell                = "Cell"
showMatType (MatComplex t)         = "Complex(" ++ showMatType t ++ ")"
showMatType (MatStruct (Known fs)) = "Struct(" ++ show (intercalate ", " fs) ++ ")"
showMatType (MatArray t (Known d)) = "[" ++ showMatType t ++ "]_" ++ show d
showMatType (MatUnknown c t)       = "Unknown(Class " ++ show c ++ ", Type " ++ show t ++ ")"

-- | Contains a dictionary for the 'Show' class.
data ShowDict a where
  ShowDict :: (Show a) => ShowDict a

-- | For a given MAT type description, compute a corresponding 'Show'
-- dictionary.
matShowDict :: MatType d a b -> ShowDict b
matShowDict MatDouble        = ShowDict
matShowDict MatInt64         = ShowDict
matShowDict MatBool          = ShowDict
matShowDict MatChar          = ShowDict
matShowDict MatCell          = ShowDict
matShowDict (MatComplex t)   = case matShowDict t of
  ShowDict -> ShowDict
matShowDict (MatStruct _)    = ShowDict
matShowDict (MatArray t _)   = case matShowDict t of
  ShowDict -> ShowDict
matShowDict (MatUnknown _ _) = ShowDict

-- | This type can be used in 'MatType' in order to represent known
-- dimensions.
newtype Known a = Known a
  deriving (Show)

-- | This type can be used in 'MatType' in order to represent unknown
-- dimensions.
data Unknown a = Unknown
  deriving (Show)

-- deriving instance (Show (d Int), Show (d [String])) => Show (MatType d a b)
-- deriving instance Eq   d => Eq   (MatType d a b)
-- deriving instance Ord  d => Ord  (MatType d a b)
instance Show (MatType Known a b) where
  show t = showMatType t

-- | Type synonym for a struct.
type MatStruct = [MatData]
  -- TODO: would [(String, MatData)] or Map String MatData be better?
  -- The type already contains the field names, but perhaps a little
  -- redundancy is ok and pragmatic here?

-- | An arbitrary 'MatType' in a wrapper.
data AnyMatType :: * where
  AnyMatType :: (Storable a) => MatType Known a b -> AnyMatType

deriving instance Show (AnyMatType)
-- deriving instance Eq   (AnyMatType)
-- deriving instance Ord  (AnyMatType)

-- | Some piece of data paired with a type annotation.
data MatData :: * where
  MatData :: MatType Known a b -> b -> MatData

-- deriving instance Show (MatData)
instance Show MatData where
  show (MatData t d) = case matShowDict t of
    ShowDict -> show d ++ " :: " ++ show t

-- TODO: Think about a proper 'Show' instance.

-- | Given a type representation, 'matParser' computes a memory parser
-- for that type.
matParser :: MatType Known a b -> MemParser a b
matParser MatDouble              = storable
matParser MatInt64               = storable
matParser MatBool                = toBool <$> storable
matParser MatChar                = castCUCharToChar <$> storable
matParser (MatComplex t)         = complex (matParser t)
matParser (MatStruct (Known d))  = do
  ptrs <- array (length d) storable
  liftIO $ forM ptrs $ matDataPtr
matParser MatCell                = do
  ptr <- storable
  liftIO (matDataPtr ptr)
matParser (MatArray t (Known n)) = array n (matParser t)
matParser (MatUnknown _ _)       = nothing

-- | Given the raw variable data structure, get the name
-- of the variable.
matVarName :: CMatVar -> IO (Maybe String)
matVarName v = maybePeek peekCString (name v)

-- | Given the raw variable data structure, get the dimensions
-- of the variable.
--
-- TODO: Can the dimensions pointer be null?
matVarDims :: CMatVar -> IO [Int]
matVarDims v = map fromIntegral <$> peekArray (fromIntegral $ rank v) (dims v)

-- | Compute the type of the represented data.
matType :: Ptr CMatVar -> CMatVar -> IO AnyMatType
matType ptr v = do
  base <- matBaseType
  ds   <- matVarDims v
  return (matComplex (matDims ds base))
  where
    matBaseType :: IO AnyMatType
    matBaseType
      | dataType v == tdouble = return $ AnyMatType MatDouble
      | dataType v == tint64  = return $ AnyMatType MatInt64
      | dataType v == tuint8 && classType v == cchar = return $ AnyMatType MatChar
      | dataType v == tuint8 && isLogical v          = return $ AnyMatType MatBool
      | dataType v == tstruct = do
          nrfields <- cMatVarGetNumberOfFields ptr
          fields   <- cMatVarGetStructFieldnames >=> peekArray (fromIntegral nrfields) >=> mapM peekCString $ ptr
          return $ AnyMatType (MatStruct (Known fields))
      | dataType v == tcell   = return $ AnyMatType MatCell
      | otherwise             = return $ AnyMatType $ MatUnknown (classType v) (dataType v)

    matDims :: [Int] -> AnyMatType -> AnyMatType
    matDims []       t = t
    matDims (d : ds) t =
      case matDims ds t of
        AnyMatType r -> AnyMatType (MatArray r (Known (fromIntegral d)))

    matComplex :: AnyMatType -> AnyMatType
    matComplex t | isComplex v =
      case t of
        AnyMatType r -> AnyMatType (MatComplex r)
                 | otherwise   = t

matTypePtr :: Ptr CMatVar -> IO AnyMatType
matTypePtr ptr = peek ptr >>= matType ptr

-- | Given a pointer to a MAT variable, provide the data in a suitable
-- Haskell structure, wrapped up with a type representation.
matData :: Ptr CMatVar -> CMatVar -> IO MatData
matData ptr v = do
  AnyMatType t <- matType ptr v
  let mp = matParser t
  x <- runMemParser mp (castPtr (theData v))
  return (MatData t x)

matDataPtr :: Ptr CMatVar -> IO MatData
matDataPtr ptr = peek ptr >>= matData ptr

data Equal :: * -> * -> * where
  Refl :: Equal a a

matTypeEqual :: MatType d a b -> MatType d' a' c -> Maybe (Equal b c)
matTypeEqual MatDouble  MatDouble = return Refl
matTypeEqual MatInt64   MatInt64  = return Refl
matTypeEqual MatBool    MatBool   = return Refl
matTypeEqual MatChar    MatChar   = return Refl
matTypeEqual MatCell    MatCell   = return Refl
matTypeEqual (MatComplex s) (MatComplex t) = do
  Refl <- matTypeEqual s t
  return Refl
matTypeEqual (MatStruct _) (MatStruct _) = return Refl
matTypeEqual (MatArray s _) (MatArray t _) = do
  Refl <- matTypeEqual s t
  return Refl
matTypeEqual _ _ = Nothing

-- | Given a value of unknown type and a type representation,
-- attempt a cast.
matCast :: MatData -> MatType d a b -> Maybe b
matCast (MatData s x) t =
  case matTypeEqual s t of
    Just Refl -> Just x
    Nothing   -> Nothing

dim2 :: MatType Unknown a b -> MatType Unknown a [[b]]
dim2 t = MatArray (MatArray t Unknown) Unknown

max_stars_repo_name: glutamate/matio

max_issues_repo_name: glutamate/matio

max_forks_repo_name: glutamate/matio


text:
import Control.Parallel

import Numeric.Container hiding (linspace)
import Numeric.LinearAlgebra.Util
import Numeric.LinearAlgebra.Data hiding (linspace)
import Numeric.LinearAlgebra.Algorithms
import Numeric.LinearAlgebra.Devel
import Data.List hiding (find)
import DSP.Basic

--For Charts
--import Graphics.Rendering.Chart hiding (Matrix, Vector)
--import Data.Colour
--import Data.Colour.Names
--import Data.Default.Class
--import Graphics.Rendering.Chart.Backend.Cairo
--import Control.Lens

--gridsize
n = 60

--index function
index :: (Int, Int) -> Int -> Double
index (i,j) n
	| i == j       =  4
	| i == j + n   = -1
	| i == j - n   = -1
	| i == j + 1   = -1
	| i == j - 1   = -1
	| otherwise    =  0

-- Laplace Matrix in 2 dimensions
t :: Matrix Double
t = buildMatrix (n*n) (n*n) (\(i,j) -> index (i,j) n)

-- Compute eigenvalues and corresponding eigenvectors
e = sort(toList(mapVector realPart (fst(eig(t)))))
psi = toColumns (mapMatrix realPart (snd(eig(t))))

--For further image processing
grid :: Vector Double
grid = fromList (linspace 0.0 1.0 n)

ygrid = fromColumns $ replicate n grid
xgrid = fromRows    $ replicate n grid
-- eigenfunction matrix of n=12 for plotting
z1grid = reshape n $ mapVector (\x -> abs(x)^2) (psi !! 12)

main :: IO ()
main = do
-- write out matrices to visualize with python
saveMatrix "xmat.txt" "%f" xgrid
saveMatrix "ymat.txt" "%f" ygrid
saveMatrix "zmat.txt" "%f" z1grid

max_stars_repo_name: davidschlegel/hartree-fock

max_issues_repo_name: davidschlegel/hartree-fock

max_forks_repo_name: davidschlegel/hartree-fock


text:
{-# LANGUAGE TemplateHaskell #-}

module Types where

import Lens.Micro.TH
import Lens.Micro
import qualified Statistics.Distribution.Normal as Stat
import qualified Statistics.Distribution as Stat

data Bandit d = Bandit d

instance (Stat.Mean d, Stat.Variance d) => Show (Bandit d) where
  show (Bandit d) = "Bandit {mean = " ++ show (Stat.mean d) ++ ", stdDev = " ++ show (Stat.stdDev d) ++ "}"

distLens :: Lens' (Bandit d) d
distLens = lens (\(Bandit d) -> d) (\_ d -> Bandit d)

meanLens :: Lens' Stat.NormalDistribution Double
meanLens = lens (\d -> Stat.mean d) (\d m -> Stat.normalDistr m (Stat.stdDev d))

-- Statistics about a given bandit.
data BanditStats = BanditStats { _totalReward :: Double
                               , _timesPulled :: Int
                               } deriving (Eq, Show)

makeLenses ''BanditStats

instance Monoid BanditStats where
  a `mappend` b = BanditStats (a ^. totalReward + b ^. totalReward) (a ^. timesPulled + b ^. timesPulled)
  mempty = BanditStats 0.0 0

type BanditDist = Stat.NormalDistribution

data SimState = SimState { _bandits :: [(BanditStats, Bandit BanditDist)] -- ^ The bandits and associated statistics
                         , _roundsCount :: Int -- ^ Number of rounds so far
                         , _regretHist :: [Double]
                         }

makeLenses ''SimState

instance Show SimState where
  show (SimState b count regretHist_) = "SimState { bandits = " ++ show b ++ ", n = " ++ show count ++ ", regret = "++ show regretHist_ ++"}"

type BanditState = (BanditStats, Bandit BanditDist)

max_stars_repo_name: beala/bandits

max_issues_repo_name: beala/bandits

max_forks_repo_name: beala/bandits


text:
{-# LANGUAGE BangPatterns, CPP, MagicHash,
             ScopedTypeVariables, UnliftedFFITypes, DeriveDataTypeable,
             DefaultSignatures, FlexibleContexts, TypeFamilies,
             MultiParamTypeClasses, CApiFFI #-}

{-# LANGUAGE Trustworthy #-}

#if __GLASGOW_HASKELL__ >= 801
{-# LANGUAGE PolyKinds #-} -- For TypeRep instances
#endif

{-# OPTIONS_GHC -fno-warn-deprecations #-}

------------------------------------------------------------------------
-- |
-- Module      :  Data.Hashable.Class
-- Copyright   :  (c) Milan Straka 2010
--                (c) Johan Tibell 2011
--                (c) Bryan O'Sullivan 2011, 2012
-- SPDX-License-Identifier : BSD-3-Clause
-- Maintainer  :  johan.tibell@gmail.com
-- Stability   :  provisional
-- Portability :  portable
--
-- This module defines a class, 'Hashable', for types that can be
-- converted to a hash value.  This class exists for the benefit of
-- hashing-based data structures.  The module provides instances for
-- most standard types.

module Data.Hashable.Class
    (
      -- * Computing hash values
      Hashable(..)
    , Hashable1(..)
    , Hashable2(..)

      -- ** Support for generics
    , genericHashWithSalt
    , genericLiftHashWithSalt
    , GHashable(..)
    , HashArgs(..)
    , Zero
    , One

      -- * Creating new instances
    , hashUsing
    , hashPtr
    , hashPtrWithSalt
    , hashByteArray
    , hashByteArrayWithSalt
    , defaultHashWithSalt
      -- * Higher Rank Functions
    , hashWithSalt1
    , hashWithSalt2
    , defaultLiftHashWithSalt
    -- * Caching hashes
    , Hashed
    , hashed
    , unhashed
    , mapHashed
    , traverseHashed
    ) where

import Control.Applicative (Const(..))
import Control.Exception (assert)
import Control.DeepSeq (NFData(rnf))
import Data.Bits (shiftL, shiftR, xor)
import qualified Data.ByteString as B
import qualified Data.ByteString.Lazy as BL
import qualified Data.ByteString.Unsafe as B
import Data.Complex (Complex(..))
import Data.Int (Int8, Int16, Int32, Int64)
import Data.List (foldl')
import Data.Ratio (Ratio, denominator, numerator)
import qualified Data.Text as T
import qualified Data.Text.Array as TA
import qualified Data.Text.Internal as T
import qualified Data.Text.Lazy as TL
import Data.Version (Version(..))
import Data.Word (Word8, Word16, Word32, Word64)
import Foreign.C (CString)
import Foreign.Marshal.Utils (with)
import Foreign.Ptr (Ptr, FunPtr, IntPtr, WordPtr, castPtr, castFunPtrToPtr, ptrToIntPtr)
import Foreign.Storable (alignment, peek, sizeOf)
import GHC.Base (ByteArray#)
import GHC.Conc (ThreadId(..))
import GHC.Prim (ThreadId#)
import System.IO.Unsafe (unsafeDupablePerformIO)
import System.Mem.StableName
import Data.Unique (Unique, hashUnique)

-- As we use qualified F.Foldable, we don't get warnings with newer base
import qualified Data.Foldable as F

#if MIN_VERSION_base(4,7,0)
import Data.Proxy (Proxy)
#endif

#if MIN_VERSION_base(4,7,0)
import Data.Fixed (Fixed(..))
#else
import Data.Fixed (Fixed)
import Unsafe.Coerce (unsafeCoerce)
#endif

#if MIN_VERSION_base(4,8,0)
import Data.Functor.Identity (Identity(..))
#endif

import GHC.Generics

#if   MIN_VERSION_base(4,10,0)
import Type.Reflection (Typeable, TypeRep, SomeTypeRep(..))
import Type.Reflection.Unsafe (typeRepFingerprint)
import GHC.Fingerprint.Type(Fingerprint(..))
#elif MIN_VERSION_base(4,8,0)
import Data.Typeable (typeRepFingerprint, Typeable, TypeRep)
import GHC.Fingerprint.Type(Fingerprint(..))
#else
import Data.Typeable.Internal (Typeable, TypeRep (..))
import GHC.Fingerprint.Type(Fingerprint(..))
#endif

import Foreign.C.Types (CInt(..))

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

#if MIN_VERSION_base(4,7,0)
import Data.Bits (finiteBitSize)
#else
import Data.Bits (bitSize)
#endif

#if !(MIN_VERSION_bytestring(0,10,0))
import qualified Data.ByteString.Lazy.Internal as BL  -- foldlChunks
#endif

#if MIN_VERSION_bytestring(0,10,4)
import qualified Data.ByteString.Short.Internal as BSI
#endif

#ifdef VERSION_ghc_bignum
import GHC.Num.BigNat (BigNat (..))
import GHC.Num.Integer (Integer (..))
import GHC.Num.Natural (Natural (..))
import GHC.Exts (Int (..), sizeofByteArray#)
#endif

#ifdef VERSION_integer_gmp

# if MIN_VERSION_integer_gmp(1,0,0)
#  define MIN_VERSION_integer_gmp_1_0_0
# endif

import GHC.Exts (Int(..))
import GHC.Integer.GMP.Internals (Integer(..))
# if defined(MIN_VERSION_integer_gmp_1_0_0)
import GHC.Exts (sizeofByteArray#)
import GHC.Integer.GMP.Internals (BigNat(BN#))
# endif
#endif

#if MIN_VERSION_base(4,8,0)
import Data.Void (Void, absurd)
import GHC.Exts (Word(..))
#ifndef VERSION_ghc_bignum
import GHC.Natural (Natural(..))
#endif
#endif

#if MIN_VERSION_base(4,9,0)
import qualified Data.List.NonEmpty as NE
import Data.Semigroup
import Data.Functor.Classes (Eq1(..),Ord1(..),Show1(..),showsUnaryWith)

import Data.Functor.Compose (Compose(..))
import qualified Data.Functor.Product as FP
import qualified Data.Functor.Sum as FS
#endif

import Data.String (IsString(..))

#if MIN_VERSION_base(4,9,0)
import Data.Kind (Type)
#else
#define Type *
#endif

#ifdef HASHABLE_RANDOM_SEED
import System.IO.Unsafe (unsafePerformIO)
#endif

#include "MachDeps.h"

infixl 0 `hashWithSalt`

------------------------------------------------------------------------
-- * Computing hash values

#ifdef HASHABLE_RANDOM_SEED
initialSeed :: Word64
initialSeed = unsafePerformIO initialSeedC
{-# NOINLINE initialSeed #-}

foreign import capi "HsHashable.h hs_hashable_init" initialSeedC :: IO Word64
#endif

-- | A default salt used in the implementation of 'hash'.
defaultSalt :: Int
#ifdef HASHABLE_RANDOM_SEED
defaultSalt = hashWithSalt defaultSalt' initialSeed
#else
defaultSalt = defaultSalt'
#endif
{-# INLINE defaultSalt #-}

defaultSalt' :: Int
#if WORD_SIZE_IN_BITS == 64
defaultSalt' = -3750763034362895579 -- 14695981039346656037 :: Int64
#else
defaultSalt' = -2128831035 -- 2166136261 :: Int32
#endif
{-# INLINE defaultSalt' #-}

-- | The class of types that can be converted to a hash value.
--
-- Minimal implementation: 'hashWithSalt'.
--
-- /Note:/ the hash is not guaranteed to be stable across
-- library versions, operating systems or architectures.
-- For stable hashing use named hashes: SHA256, CRC32 etc.
--
-- If you are looking for 'Hashable' instance in @time@ package,
-- check [time-compat](https://hackage.haskell.org/package/time-compat)
--
class Hashable a where
    -- | Return a hash value for the argument, using the given salt.
    --
    -- The general contract of 'hashWithSalt' is:
    --
    --  * If two values are equal according to the '==' method, then
    --    applying the 'hashWithSalt' method on each of the two values
    --    /must/ produce the same integer result if the same salt is
    --    used in each case.
    --
    --  * It is /not/ required that if two values are unequal
    --    according to the '==' method, then applying the
    --    'hashWithSalt' method on each of the two values must produce
    --    distinct integer results. However, the programmer should be
    --    aware that producing distinct integer results for unequal
    --    values may improve the performance of hashing-based data
    --    structures.
    --
    --  * This method can be used to compute different hash values for
    --    the same input by providing a different salt in each
    --    application of the method. This implies that any instance
    --    that defines 'hashWithSalt' /must/ make use of the salt in
    --    its implementation.
    --
    --  * 'hashWithSalt' may return negative 'Int' values.
    --
    hashWithSalt :: Int -> a -> Int

    -- | Like 'hashWithSalt', but no salt is used. The default
    -- implementation uses 'hashWithSalt' with some default salt.
    -- Instances might want to implement this method to provide a more
    -- efficient implementation than the default implementation.
    hash :: a -> Int
    hash = hashWithSalt defaultSalt

    default hashWithSalt :: (Generic a, GHashable Zero (Rep a)) => Int -> a -> Int
    hashWithSalt = genericHashWithSalt
    {-# INLINE hashWithSalt #-}

-- | Generic 'hashWithSalt'.
--
-- @since 1.3.0.0
genericHashWithSalt :: (Generic a, GHashable Zero (Rep a)) => Int -> a -> Int
genericHashWithSalt = \salt -> ghashWithSalt HashArgs0 salt . from
{-# INLINE genericHashWithSalt #-}

data Zero
data One

data family HashArgs arity a :: Type
data instance HashArgs Zero a = HashArgs0
newtype instance HashArgs One  a = HashArgs1 (Int -> a -> Int)

-- | The class of types that can be generically hashed.
class GHashable arity f where
    ghashWithSalt :: HashArgs arity a -> Int -> f a -> Int

class Hashable1 t where
    -- | Lift a hashing function through the type constructor.
    liftHashWithSalt :: (Int -> a -> Int) -> Int -> t a -> Int

    default liftHashWithSalt :: (Generic1 t, GHashable One (Rep1 t)) => (Int -> a -> Int) -> Int -> t a -> Int
    liftHashWithSalt = genericLiftHashWithSalt
    {-# INLINE liftHashWithSalt #-}

-- | Generic 'liftHashWithSalt'.
--
-- @since 1.3.0.0
genericLiftHashWithSalt :: (Generic1 t, GHashable One (Rep1 t)) => (Int -> a -> Int) -> Int -> t a -> Int
genericLiftHashWithSalt = \h salt -> ghashWithSalt (HashArgs1 h) salt . from1
{-# INLINE genericLiftHashWithSalt #-}

class Hashable2 t where
    -- | Lift a hashing function through the binary type constructor.
    liftHashWithSalt2 :: (Int -> a -> Int) -> (Int -> b -> Int) -> Int -> t a b -> Int

-- | Lift the 'hashWithSalt' function through the type constructor.
--
-- > hashWithSalt1 = liftHashWithSalt hashWithSalt
hashWithSalt1 :: (Hashable1 f, Hashable a) => Int -> f a -> Int
hashWithSalt1 = liftHashWithSalt hashWithSalt

-- | Lift the 'hashWithSalt' function through the type constructor.
--
-- > hashWithSalt2 = liftHashWithSalt2 hashWithSalt hashWithSalt
hashWithSalt2 :: (Hashable2 f, Hashable a, Hashable b) => Int -> f a b -> Int
hashWithSalt2 = liftHashWithSalt2 hashWithSalt hashWithSalt

-- | Lift the 'hashWithSalt' function halfway through the type constructor.
-- This function makes a suitable default implementation of 'liftHashWithSalt',
-- given that the type constructor @t@ in question can unify with @f a@.
defaultLiftHashWithSalt :: (Hashable2 f, Hashable a) => (Int -> b -> Int) -> Int -> f a b -> Int
defaultLiftHashWithSalt h = liftHashWithSalt2 hashWithSalt h

-- | Since we support a generic implementation of 'hashWithSalt' we
-- cannot also provide a default implementation for that method for
-- the non-generic instance use case. Instead we provide
-- 'defaultHashWith'.
defaultHashWithSalt :: Hashable a => Int -> a -> Int
defaultHashWithSalt salt x = salt `combine` hash x

-- | Transform a value into a 'Hashable' value, then hash the
-- transformed value using the given salt.
--
-- This is a useful shorthand in cases where a type can easily be
-- mapped to another type that is already an instance of 'Hashable'.
-- Example:
--
-- > data Foo = Foo | Bar
-- >          deriving (Enum)
-- >
-- > instance Hashable Foo where
-- >     hashWithSalt = hashUsing fromEnum
--
-- @since 1.2.0.0
hashUsing :: (Hashable b) =>
             (a -> b)           -- ^ Transformation function.
          -> Int                -- ^ Salt.
          -> a                  -- ^ Value to transform.
          -> Int
hashUsing f salt x = hashWithSalt salt (f x)
{-# INLINE hashUsing #-}

instance Hashable Int where
    hash = id
    hashWithSalt = defaultHashWithSalt

instance Hashable Int8 where
    hash = fromIntegral
    hashWithSalt = defaultHashWithSalt

instance Hashable Int16 where
    hash = fromIntegral
    hashWithSalt = defaultHashWithSalt

instance Hashable Int32 where
    hash = fromIntegral
    hashWithSalt = defaultHashWithSalt

instance Hashable Int64 where
    hash n
#if MIN_VERSION_base(4,7,0)
        | finiteBitSize (undefined :: Int) == 64 = fromIntegral n
#else
        | bitSize (undefined :: Int) == 64 = fromIntegral n
#endif
        | otherwise = fromIntegral (fromIntegral n `xor`
                                   (fromIntegral n `shiftR` 32 :: Word64))
    hashWithSalt = defaultHashWithSalt

instance Hashable Word where
    hash = fromIntegral
    hashWithSalt = defaultHashWithSalt

instance Hashable Word8 where
    hash = fromIntegral
    hashWithSalt = defaultHashWithSalt

instance Hashable Word16 where
    hash = fromIntegral
    hashWithSalt = defaultHashWithSalt

instance Hashable Word32 where
    hash = fromIntegral
    hashWithSalt = defaultHashWithSalt

instance Hashable Word64 where
    hash n
#if MIN_VERSION_base(4,7,0)
        | finiteBitSize (undefined :: Int) == 64 = fromIntegral n
#else
        | bitSize (undefined :: Int) == 64 = fromIntegral n
#endif
        | otherwise = fromIntegral (n `xor` (n `shiftR` 32))
    hashWithSalt = defaultHashWithSalt

instance Hashable () where
    hash = fromEnum
    hashWithSalt = defaultHashWithSalt

instance Hashable Bool where
    hash = fromEnum
    hashWithSalt = defaultHashWithSalt

instance Hashable Ordering where
    hash = fromEnum
    hashWithSalt = defaultHashWithSalt

instance Hashable Char where
    hash = fromEnum
    hashWithSalt = defaultHashWithSalt

#if defined(MIN_VERSION_integer_gmp_1_0_0) || defined(VERSION_ghc_bignum)
instance Hashable BigNat where
    hashWithSalt salt (BN# ba) = hashByteArrayWithSalt ba 0 numBytes salt
                                 `hashWithSalt` size
      where
        size     = numBytes `quot` SIZEOF_HSWORD
        numBytes = I# (sizeofByteArray# ba)
#endif

#if MIN_VERSION_base(4,8,0)
instance Hashable Natural where
# if defined(VERSION_ghc_bignum)
    hash (NS n)   = hash (W# n)
    hash (NB bn)  = hash (BN# bn)

    hashWithSalt salt (NS n)  = hashWithSalt salt (W# n)
    hashWithSalt salt (NB bn) = hashWithSalt salt (BN# bn)
# else
# if defined(MIN_VERSION_integer_gmp_1_0_0)
    hash (NatS# n)   = hash (W# n)
    hash (NatJ# bn)  = hash bn

    hashWithSalt salt (NatS# n)   = hashWithSalt salt (W# n)
    hashWithSalt salt (NatJ# bn)  = hashWithSalt salt bn
# else
    hash (Natural n) = hash n

    hashWithSalt salt (Natural n) = hashWithSalt salt n
# endif
# endif
#endif

instance Hashable Integer where
#if defined(VERSION_ghc_bignum)
    hash (IS n)  = I# n
    hash (IP bn) = hash (BN# bn)
    hash (IN bn) = negate (hash (BN# bn))

    hashWithSalt salt (IS n)  = hashWithSalt salt (I# n)
    hashWithSalt salt (IP bn) = hashWithSalt salt (BN# bn)
    hashWithSalt salt (IN bn) = negate (hashWithSalt salt (BN# bn))
#else
#if defined(VERSION_integer_gmp)
# if defined(MIN_VERSION_integer_gmp_1_0_0)
    hash (S# n)   = (I# n)
    hash (Jp# bn) = hash bn
    hash (Jn# bn) = negate (hash bn)

    hashWithSalt salt (S# n)   = hashWithSalt salt (I# n)
    hashWithSalt salt (Jp# bn) = hashWithSalt salt bn
    hashWithSalt salt (Jn# bn) = negate (hashWithSalt salt bn)
# else
    hash (S# int) = I# int
    hash n@(J# size# byteArray)
        | n >= minInt && n <= maxInt = fromInteger n :: Int
        | otherwise = let size = I# size#
                          numBytes = SIZEOF_HSWORD * abs size
                      in hashByteArrayWithSalt byteArray 0 numBytes defaultSalt
                         `hashWithSalt` size
      where minInt = fromIntegral (minBound :: Int)
            maxInt = fromIntegral (maxBound :: Int)

    hashWithSalt salt (S# n) = hashWithSalt salt (I# n)
    hashWithSalt salt n@(J# size# byteArray)
        | n >= minInt && n <= maxInt = hashWithSalt salt (fromInteger n :: Int)
        | otherwise = let size = I# size#
                          numBytes = SIZEOF_HSWORD * abs size
                      in hashByteArrayWithSalt byteArray 0 numBytes salt
                         `hashWithSalt` size
      where minInt = fromIntegral (minBound :: Int)
            maxInt = fromIntegral (maxBound :: Int)
# endif
#else
    hashWithSalt salt = foldl' hashWithSalt salt . go
      where
        go n | inBounds n = [fromIntegral n :: Int]
             | otherwise   = fromIntegral n : go (n `shiftR` WORD_SIZE_IN_BITS)
        maxInt = fromIntegral (maxBound :: Int)
        inBounds x = x >= fromIntegral (minBound :: Int) && x <= maxInt
#endif
#endif

instance Hashable a => Hashable (Complex a) where
    {-# SPECIALIZE instance Hashable (Complex Double) #-}
    {-# SPECIALIZE instance Hashable (Complex Float)  #-}
    hash (r :+ i) = hash r `hashWithSalt` i
    hashWithSalt = hashWithSalt1
instance Hashable1 Complex where
    liftHashWithSalt h s (r :+ i) = s `h` r `h` i

#if MIN_VERSION_base(4,9,0)
-- Starting with base-4.9, numerator/denominator don't need 'Integral' anymore
instance Hashable a => Hashable (Ratio a) where
#else
instance (Integral a, Hashable a) => Hashable (Ratio a) where
#endif
    {-# SPECIALIZE instance Hashable (Ratio Integer) #-}
    hash a = hash (numerator a) `hashWithSalt` denominator a
    hashWithSalt s a = s `hashWithSalt` numerator a `hashWithSalt` denominator a

-- | __Note__: prior to @hashable-1.3.0.0@, @hash 0.0 /= hash (-0.0)@
--
-- The 'hash' of NaN is not well defined.
--
-- @since 1.3.0.0
instance Hashable Float where
    hash x
        | x == -0.0 || x == 0.0 = 0 -- see note in 'Hashable Double'
        | isIEEE x =
            assert (sizeOf x >= sizeOf (0::Word32) &&
                    alignment x >= alignment (0::Word32)) $
            hash ((unsafeDupablePerformIO $ with x $ peek . castPtr) :: Word32)
        | otherwise = hash (show x)
    hashWithSalt = defaultHashWithSalt

-- | __Note__: prior to @hashable-1.3.0.0@, @hash 0.0 /= hash (-0.0)@
--
-- The 'hash' of NaN is not well defined.
--
-- @since 1.3.0.0
instance Hashable Double where
    hash x
        | x == -0.0 || x == 0.0 = 0 -- s.t. @hash -0.0 == hash 0.0@ ; see #173
        | isIEEE x =
            assert (sizeOf x >= sizeOf (0::Word64) &&
                    alignment x >= alignment (0::Word64)) $
            hash ((unsafeDupablePerformIO $ with x $ peek . castPtr) :: Word64)
        | otherwise = hash (show x)
    hashWithSalt = defaultHashWithSalt

-- | A value with bit pattern (01)* (or 5* in hexa), for any size of Int.
-- It is used as data constructor distinguisher. GHC computes its value during
-- compilation.
distinguisher :: Int
distinguisher = fromIntegral $ (maxBound :: Word) `quot` 3
{-# INLINE distinguisher #-}

instance Hashable a => Hashable (Maybe a) where
    hash Nothing = 0
    hash (Just a) = distinguisher `hashWithSalt` a
    hashWithSalt = hashWithSalt1

instance Hashable1 Maybe where
    liftHashWithSalt _ s Nothing = s `combine` 0
    liftHashWithSalt h s (Just a) = s `combine` distinguisher `h` a

instance (Hashable a, Hashable b) => Hashable (Either a b) where
    hash (Left a)  = 0 `hashWithSalt` a
    hash (Right b) = distinguisher `hashWithSalt` b
    hashWithSalt = hashWithSalt1

instance Hashable a => Hashable1 (Either a) where
    liftHashWithSalt = defaultLiftHashWithSalt

instance Hashable2 Either where
    liftHashWithSalt2 h _ s (Left a) = s `combine` 0 `h` a
    liftHashWithSalt2 _ h s (Right b) = s `combine` distinguisher `h` b

instance (Hashable a1, Hashable a2) => Hashable (a1, a2) where
    hash (a1, a2) = hash a1 `hashWithSalt` a2
    hashWithSalt = hashWithSalt1

instance Hashable a1 => Hashable1 ((,) a1) where
    liftHashWithSalt = defaultLiftHashWithSalt

instance Hashable2 (,) where
    liftHashWithSalt2 h1 h2 s (a1, a2) = s `h1` a1 `h2` a2

instance (Hashable a1, Hashable a2, Hashable a3) => Hashable (a1, a2, a3) where
    hash (a1, a2, a3) = hash a1 `hashWithSalt` a2 `hashWithSalt` a3
    hashWithSalt = hashWithSalt1

instance (Hashable a1, Hashable a2) => Hashable1 ((,,) a1 a2) where
    liftHashWithSalt = defaultLiftHashWithSalt

instance Hashable a1 => Hashable2 ((,,) a1) where
    liftHashWithSalt2 h1 h2 s (a1, a2, a3) =
      (s `hashWithSalt` a1) `h1` a2 `h2` a3

instance (Hashable a1, Hashable a2, Hashable a3, Hashable a4) =>
         Hashable (a1, a2, a3, a4) where
    hash (a1, a2, a3, a4) = hash a1 `hashWithSalt` a2
                            `hashWithSalt` a3 `hashWithSalt` a4
    hashWithSalt = hashWithSalt1

instance (Hashable a1, Hashable a2, Hashable a3) => Hashable1 ((,,,) a1 a2 a3) where
    liftHashWithSalt = defaultLiftHashWithSalt

instance (Hashable a1, Hashable a2) => Hashable2 ((,,,) a1 a2) where
    liftHashWithSalt2 h1 h2 s (a1, a2, a3, a4) =
      (s `hashWithSalt` a1 `hashWithSalt` a2) `h1` a3 `h2` a4

instance (Hashable a1, Hashable a2, Hashable a3, Hashable a4, Hashable a5)
      => Hashable (a1, a2, a3, a4, a5) where
    hash (a1, a2, a3, a4, a5) =
        hash a1 `hashWithSalt` a2 `hashWithSalt` a3
        `hashWithSalt` a4 `hashWithSalt` a5
    hashWithSalt = hashWithSalt1

instance (Hashable a1, Hashable a2, Hashable a3,
          Hashable a4) => Hashable1 ((,,,,) a1 a2 a3 a4) where
    liftHashWithSalt = defaultLiftHashWithSalt

instance (Hashable a1, Hashable a2, Hashable a3)
      => Hashable2 ((,,,,) a1 a2 a3) where
    liftHashWithSalt2 h1 h2 s (a1, a2, a3, a4, a5) =
      (s `hashWithSalt` a1 `hashWithSalt` a2
         `hashWithSalt` a3) `h1` a4 `h2` a5

instance (Hashable a1, Hashable a2, Hashable a3, Hashable a4, Hashable a5,
          Hashable a6) => Hashable (a1, a2, a3, a4, a5, a6) where
    hash (a1, a2, a3, a4, a5, a6) =
        hash a1 `hashWithSalt` a2 `hashWithSalt` a3
        `hashWithSalt` a4 `hashWithSalt` a5 `hashWithSalt` a6
    hashWithSalt = hashWithSalt1

instance (Hashable a1, Hashable a2, Hashable a3, Hashable a4,
          Hashable a5) => Hashable1 ((,,,,,) a1 a2 a3 a4 a5) where
    liftHashWithSalt = defaultLiftHashWithSalt

instance (Hashable a1, Hashable a2, Hashable a3,
          Hashable a4) => Hashable2 ((,,,,,) a1 a2 a3 a4) where
    liftHashWithSalt2 h1 h2 s (a1, a2, a3, a4, a5, a6) =
      (s `hashWithSalt` a1 `hashWithSalt` a2 `hashWithSalt` a3
         `hashWithSalt` a4) `h1` a5 `h2` a6

instance (Hashable a1, Hashable a2, Hashable a3, Hashable a4, Hashable a5,
          Hashable a6, Hashable a7) =>
         Hashable (a1, a2, a3, a4, a5, a6, a7) where
    hash (a1, a2, a3, a4, a5, a6, a7) =
        hash a1 `hashWithSalt` a2 `hashWithSalt` a3
        `hashWithSalt` a4 `hashWithSalt` a5 `hashWithSalt` a6 `hashWithSalt` a7
    hashWithSalt s (a1, a2, a3, a4, a5, a6, a7) =
        s `hashWithSalt` a1 `hashWithSalt` a2 `hashWithSalt` a3
        `hashWithSalt` a4 `hashWithSalt` a5 `hashWithSalt` a6 `hashWithSalt` a7

instance (Hashable a1, Hashable a2, Hashable a3, Hashable a4, Hashable a5, Hashable a6) => Hashable1 ((,,,,,,) a1 a2 a3 a4 a5 a6) where
    liftHashWithSalt = defaultLiftHashWithSalt

instance (Hashable a1, Hashable a2, Hashable a3, Hashable a4,
          Hashable a5) => Hashable2 ((,,,,,,) a1 a2 a3 a4 a5) where
    liftHashWithSalt2 h1 h2 s (a1, a2, a3, a4, a5, a6, a7) =
      (s `hashWithSalt` a1 `hashWithSalt` a2 `hashWithSalt` a3
         `hashWithSalt` a4 `hashWithSalt` a5) `h1` a6 `h2` a7

instance Hashable (StableName a) where
    hash = hashStableName
    hashWithSalt = defaultHashWithSalt

-- Auxillary type for Hashable [a] definition
data SPInt = SP !Int !Int

instance Hashable a => Hashable [a] where
    {-# SPECIALIZE instance Hashable [Char] #-}
    hashWithSalt = hashWithSalt1

instance Hashable1 [] where
    liftHashWithSalt h salt arr = finalise (foldl' step (SP salt 0) arr)
      where
        finalise (SP s l) = hashWithSalt s l
        step (SP s l) x   = SP (h s x) (l + 1)

instance Hashable B.ByteString where
    hashWithSalt salt bs = unsafeDupablePerformIO $
                           B.unsafeUseAsCStringLen bs $ \(p, len) ->
                           hashPtrWithSalt p (fromIntegral len) salt

instance Hashable BL.ByteString where
    hashWithSalt = BL.foldlChunks hashWithSalt

#if MIN_VERSION_bytestring(0,10,4)
instance Hashable BSI.ShortByteString where
    hashWithSalt salt sbs@(BSI.SBS ba) =
        hashByteArrayWithSalt ba 0 (BSI.length sbs) salt
#endif

instance Hashable T.Text where
    hashWithSalt salt (T.Text arr off len) =
        hashByteArrayWithSalt (TA.aBA arr) (off `shiftL` 1) (len `shiftL` 1)
        salt

instance Hashable TL.Text where
    hashWithSalt = TL.foldlChunks hashWithSalt

-- | Compute the hash of a ThreadId.
hashThreadId :: ThreadId -> Int
hashThreadId (ThreadId t) = hash (fromIntegral (getThreadId t) :: Int)

-- this cannot be capi, as GHC panics.
foreign import ccall unsafe "rts_getThreadId" getThreadId
    :: ThreadId# -> CInt

instance Hashable ThreadId where
    hash = hashThreadId
    hashWithSalt = defaultHashWithSalt

instance Hashable (Ptr a) where
    hashWithSalt salt p = hashWithSalt salt $ ptrToIntPtr p

instance Hashable (FunPtr a) where
    hashWithSalt salt p = hashWithSalt salt $ castFunPtrToPtr p

instance Hashable IntPtr where
    hash n = fromIntegral n
    hashWithSalt = defaultHashWithSalt

instance Hashable WordPtr where
    hash n = fromIntegral n
    hashWithSalt = defaultHashWithSalt

----------------------------------------------------------------------------
-- Fingerprint & TypeRep instances

-- | @since 1.3.0.0
instance Hashable Fingerprint where
    hash (Fingerprint x _) = fromIntegral x
    hashWithSalt = defaultHashWithSalt
    {-# INLINE hash #-}

#if MIN_VERSION_base(4,10,0)

hashTypeRep :: Type.Reflection.TypeRep a -> Int
hashTypeRep tr =
    let Fingerprint x _ = typeRepFingerprint tr in fromIntegral x

instance Hashable Type.Reflection.SomeTypeRep where
    hash (Type.Reflection.SomeTypeRep r) = hashTypeRep r
    hashWithSalt = defaultHashWithSalt
    {-# INLINE hash #-}

instance Hashable (Type.Reflection.TypeRep a) where
    hash = hashTypeRep
    hashWithSalt = defaultHashWithSalt
    {-# INLINE hash #-}

#else

-- | Compute the hash of a TypeRep, in various GHC versions we can do this quickly.
hashTypeRep :: TypeRep -> Int
{-# INLINE hashTypeRep #-}
#if   MIN_VERSION_base(4,8,0)
-- Fingerprint is just the MD5, so taking any Int from it is fine
hashTypeRep tr = let Fingerprint x _ = typeRepFingerprint tr in fromIntegral x
#else
-- Fingerprint is just the MD5, so taking any Int from it is fine
hashTypeRep (TypeRep (Fingerprint x _) _ _) = fromIntegral x
#endif

instance Hashable TypeRep where
    hash = hashTypeRep
    hashWithSalt = defaultHashWithSalt
    {-# INLINE hash #-}

#endif

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

#if MIN_VERSION_base(4,8,0)
instance Hashable Void where
    hashWithSalt _ = absurd
#endif

-- | Compute a hash value for the content of this pointer.
hashPtr :: Ptr a      -- ^ pointer to the data to hash
        -> Int        -- ^ length, in bytes
        -> IO Int     -- ^ hash value
hashPtr p len = hashPtrWithSalt p len defaultSalt

-- | Compute a hash value for the content of this pointer, using an
-- initial salt.
--
-- This function can for example be used to hash non-contiguous
-- segments of memory as if they were one contiguous segment, by using
-- the output of one hash as the salt for the next.
hashPtrWithSalt :: Ptr a   -- ^ pointer to the data to hash
                -> Int     -- ^ length, in bytes
                -> Int     -- ^ salt
                -> IO Int  -- ^ hash value
hashPtrWithSalt p len salt =
    fromIntegral `fmap` c_hashCString (castPtr p) (fromIntegral len)
    (fromIntegral salt)

foreign import capi unsafe "HsHashable.h hashable_fnv_hash" c_hashCString
#if WORD_SIZE_IN_BITS == 64
    :: CString -> Int64 -> Int64 -> IO Word64
#else
    :: CString -> Int32 -> Int32 -> IO Word32
#endif

-- | Compute a hash value for the content of this 'ByteArray#',
-- beginning at the specified offset, using specified number of bytes.
hashByteArray :: ByteArray#  -- ^ data to hash
              -> Int         -- ^ offset, in bytes
              -> Int         -- ^ length, in bytes
              -> Int         -- ^ hash value
hashByteArray ba0 off len = hashByteArrayWithSalt ba0 off len defaultSalt
{-# INLINE hashByteArray #-}

-- | Compute a hash value for the content of this 'ByteArray#', using
-- an initial salt.
--
-- This function can for example be used to hash non-contiguous
-- segments of memory as if they were one contiguous segment, by using
-- the output of one hash as the salt for the next.
hashByteArrayWithSalt
    :: ByteArray#  -- ^ data to hash
    -> Int         -- ^ offset, in bytes
    -> Int         -- ^ length, in bytes
    -> Int         -- ^ salt
    -> Int         -- ^ hash value
hashByteArrayWithSalt ba !off !len !h =
    fromIntegral $ c_hashByteArray ba (fromIntegral off) (fromIntegral len)
    (fromIntegral h)

#if __GLASGOW_HASKELL__ >= 802
foreign import capi unsafe "HsHashable.h hashable_fnv_hash_offset" c_hashByteArray
#else
foreign import ccall unsafe "hashable_fnv_hash_offset" c_hashByteArray
#endif
#if WORD_SIZE_IN_BITS == 64
    :: ByteArray# -> Int64 -> Int64 -> Int64 -> Word64
#else
    :: ByteArray# -> Int32 -> Int32 -> Int32 -> Word32
#endif

-- | Combine two given hash values.  'combine' has zero as a left
-- identity.
combine :: Int -> Int -> Int
combine h1 h2 = (h1 * 16777619) `xor` h2

instance Hashable Unique where
    hash = hashUnique
    hashWithSalt = defaultHashWithSalt

instance Hashable Version where
    hashWithSalt salt (Version branch tags) =
        salt `hashWithSalt` branch `hashWithSalt` tags

#if MIN_VERSION_base(4,7,0)
instance Hashable (Fixed a) where
    hashWithSalt salt (MkFixed i) = hashWithSalt salt i
-- Using hashWithSalt1 would cause needless constraint
instance Hashable1 Fixed where
    liftHashWithSalt _ salt (MkFixed i) = hashWithSalt salt i
#else
instance Hashable (Fixed a) where
    hashWithSalt salt x = hashWithSalt salt (unsafeCoerce x :: Integer)
#endif

#if MIN_VERSION_base(4,8,0)
instance Hashable a => Hashable (Identity a) where
    hashWithSalt = hashWithSalt1
instance Hashable1 Identity where
    liftHashWithSalt h salt (Identity x) = h salt x
#endif

-- Using hashWithSalt1 would cause needless constraint
instance Hashable a => Hashable (Const a b) where
    hashWithSalt salt (Const x) = hashWithSalt salt x

instance Hashable a => Hashable1 (Const a) where
    liftHashWithSalt = defaultLiftHashWithSalt

instance Hashable2 Const where
    liftHashWithSalt2 f _ salt (Const x) = f salt x

#if MIN_VERSION_base(4,7,0)
instance Hashable (Proxy a) where
    hash _ = 0
    hashWithSalt s _ = s

instance Hashable1 Proxy where
    liftHashWithSalt _ s _ = s
#endif

-- instances formerly provided by 'semigroups' package
#if MIN_VERSION_base(4,9,0)
instance Hashable a => Hashable (NE.NonEmpty a) where
    hashWithSalt p (a NE.:| as) = p `hashWithSalt` a `hashWithSalt` as

-- | @since 1.3.1.0
instance Hashable1 NE.NonEmpty where
    liftHashWithSalt h salt (a NE.:| as) = liftHashWithSalt h (h salt a) as

instance Hashable a => Hashable (Min a) where
    hashWithSalt p (Min a) = hashWithSalt p a

-- | @since 1.3.1.0
instance Hashable1 Min where liftHashWithSalt h salt (Min a) = h salt a

instance Hashable a => Hashable (Max a) where
    hashWithSalt p (Max a) = hashWithSalt p a

-- | @since 1.3.1.0
instance Hashable1 Max where liftHashWithSalt h salt (Max a) = h salt a

-- | __Note__: Prior to @hashable-1.3.0.0@ the hash computation included the second argument of 'Arg' which wasn't consistent with its 'Eq' instance.
--
-- @since 1.3.0.0
instance Hashable a => Hashable (Arg a b) where
    hashWithSalt p (Arg a _) = hashWithSalt p a

instance Hashable a => Hashable (First a) where
    hashWithSalt p (First a) = hashWithSalt p a

-- | @since 1.3.1.0
instance Hashable1 First where liftHashWithSalt h salt (First a) = h salt a

instance Hashable a => Hashable (Last a) where
    hashWithSalt p (Last a) = hashWithSalt p a

-- | @since 1.3.1.0
instance Hashable1 Last where liftHashWithSalt h salt (Last a) = h salt a

instance Hashable a => Hashable (WrappedMonoid a) where
    hashWithSalt p (WrapMonoid a) = hashWithSalt p a

-- | @since 1.3.1.0
instance Hashable1 WrappedMonoid where liftHashWithSalt h salt (WrapMonoid a) = h salt a

#if !MIN_VERSION_base(4,16,0)
instance Hashable a => Hashable (Option a) where
    hashWithSalt p (Option a) = hashWithSalt p a

-- | @since 1.3.1.0
instance Hashable1 Option where liftHashWithSalt h salt (Option a) = liftHashWithSalt h salt a
#endif
#endif

-- instances for @Data.Functor.{Product,Sum,Compose}@, present
-- in base-4.9 and onward.
#if MIN_VERSION_base(4,9,0)
-- | In general, @hash (Compose x) ≠ hash x@. However, @hashWithSalt@ satisfies
-- its variant of this equivalence.
instance (Hashable1 f, Hashable1 g, Hashable a) => Hashable (Compose f g a) where
    hashWithSalt = hashWithSalt1

instance (Hashable1 f, Hashable1 g) => Hashable1 (Compose f g) where
    liftHashWithSalt h s = liftHashWithSalt (liftHashWithSalt h) s . getCompose

instance (Hashable1 f, Hashable1 g) => Hashable1 (FP.Product f g) where
    liftHashWithSalt h s (FP.Pair a b) = liftHashWithSalt h (liftHashWithSalt h s a) b

instance (Hashable1 f, Hashable1 g, Hashable a) => Hashable (FP.Product f g a) where
    hashWithSalt = hashWithSalt1

instance (Hashable1 f, Hashable1 g) => Hashable1 (FS.Sum f g) where
    liftHashWithSalt h s (FS.InL a) = liftHashWithSalt h (s `combine` 0) a
    liftHashWithSalt h s (FS.InR a) = liftHashWithSalt h (s `combine` distinguisher) a

instance (Hashable1 f, Hashable1 g, Hashable a) => Hashable (FS.Sum f g a) where
    hashWithSalt = hashWithSalt1
#endif

-- | A hashable value along with the result of the 'hash' function.
data Hashed a = Hashed a {-# UNPACK #-} !Int
  deriving (Typeable)

-- | Wrap a hashable value, caching the 'hash' function result.
hashed :: Hashable a => a -> Hashed a
hashed a = Hashed a (hash a)

-- | Unwrap hashed value.
unhashed :: Hashed a -> a
unhashed (Hashed a _) = a

-- | Uses precomputed hash to detect inequality faster
instance Eq a => Eq (Hashed a) where
  Hashed a ha == Hashed b hb = ha == hb && a == b

instance Ord a => Ord (Hashed a) where
  Hashed a _ `compare` Hashed b _ = a `compare` b

instance Show a => Show (Hashed a) where
  showsPrec d (Hashed a _) = showParen (d > 10) $
    showString "hashed" . showChar ' ' . showsPrec 11 a

instance Hashable (Hashed a) where
  hashWithSalt = defaultHashWithSalt
  hash (Hashed _ h) = h

-- This instance is a little unsettling. It is unusal for
-- 'liftHashWithSalt' to ignore its first argument when a
-- value is actually available for it to work on.
instance Hashable1 Hashed where
  liftHashWithSalt _ s (Hashed _ h) = defaultHashWithSalt s h

instance (IsString a, Hashable a) => IsString (Hashed a) where
  fromString s = let r = fromString s in Hashed r (hash r)

instance F.Foldable Hashed where
  foldr f acc (Hashed a _) = f a acc

instance NFData a => NFData (Hashed a) where
  rnf = rnf . unhashed

-- | 'Hashed' cannot be 'Functor'
mapHashed :: Hashable b => (a -> b) -> Hashed a -> Hashed b
mapHashed f (Hashed a _) = hashed (f a)

-- | 'Hashed' cannot be 'Traversable'
traverseHashed :: (Hashable b, Functor f) => (a -> f b) -> Hashed a -> f (Hashed b)
traverseHashed f (Hashed a _) = fmap hashed (f a)

-- instances for @Data.Functor.Classes@ higher rank typeclasses
-- in base-4.9 and onward.
#if MIN_VERSION_base(4,9,0)
instance Eq1 Hashed where
  liftEq f (Hashed a ha) (Hashed b hb) = ha == hb && f a b

instance Ord1 Hashed where
  liftCompare f (Hashed a _) (Hashed b _) = f a b

instance Show1 Hashed where
  liftShowsPrec sp _ d (Hashed a _) = showsUnaryWith sp "hashed" d a
#endif

max_stars_repo_name: SupercedeTech/hashable

max_issues_repo_name: SupercedeTech/hashable

max_forks_repo_name: SupercedeTech/hashable


text:
-- PCA.hs
-- Compiles to library.
-- This library contains the functions needed to perform principle
-- component anlysis on vectorized data.
module PCA
( Hyperplane (..)
, mean
, saveHyperplane
, loadHyperplane
, linRegression
, distance
) where

import Numeric.LinearAlgebra.HMatrix
import Numeric.LinearAlgebra.Data
import Data.Packed.Matrix
import Data.Packed.Vector

data Hyperplane = Hyperplane Int (Vector Double) (Matrix Double) deriving (Show, Read)
-- First argument is the dimension of the ambient space.
-- Second argument is some vector in the hyperplane.
-- Rows of third argument form the basis of the parallel
-- linear subspace.
-- The Hyperplane is the coset of the second arg mod the third arg.

saveHyperplane :: FilePath -> Hyperplane -> IO ()
-- Saves a hyperplane in machine-readable plaintext format.
saveHyperplane f h = writeFile f (show h)

loadHyperplane :: FilePath -> IO Hyperplane
-- Loads a hyperplane that was saved in the syntax of `save`.
loadHyperplane f = readFile f >>= (return . read)

mean :: Matrix Double -> Vector Double
-- Returns the average of the rows of a matrix.
mean x = scale n . sum . toRows $ x
  where
  n = 1 / (fromIntegral . rows $ x) :: Double

normalize :: Matrix Double -> Matrix Double
-- Subtracts the average of/from the rows of a matrix.
normalize x = x - (fromRows . replicate (rows x) . mean $ x)

linRegression :: Int -> Matrix Double -> Hyperplane
-- Finds the best-fit hyperplane of the rows of a matrix.
-- First argument is the threshold to use a singular value.
-- Second argument is your data, rows are individual datapoints.
-- The first argument controls the size of the hyperplane.
-- If the first argument is too small, the hyperplane will just be
-- the whole space. If the first argument is too big, the hyper-
-- plane will just be a single point.
linRegression numsv inmat = Hyperplane n m rows
  where
  n = cols inmat
  m = mean inmat
  normalizedmat = normalize inmat
  (_,rows) = trimSVDRight numsv normalizedmat

trimSVDRight :: Int -> Matrix Double -> (Vector Double, Matrix Double)
-- Returns partial diagonal and right-side factor from SVD of
-- a matrix.
trimSVDRight numsv inmat = (v, outmat)
  where
  n' = min numsv (dim v1 - 1)
  (v1,m1) = rightSV inmat
  v = fromList . take n' . toList $ v1
  outmat = trans $ takeColumns (dim v) m1

distance :: Hyperplane -> Vector Double -> Double
-- Finds the distance from a vector to a hyperplane.
distance (Hyperplane n hv m) v
  | cols m == 0 = norm_2 v'
  | rows m == 0 = norm_2 v'
  | otherwise   = norm_2 $ distancevec
  where
  v' = v - hv
  distancevec = let a = m #> v'
                    b = trans m #> a in
                a `seq` b `seq` (b - v')


text:
{-|
Module      : Data.Meta.Multiarm
Description : Reweight multi-arm studies
Copyright   : (c) Thodoris Papakonstantinou, 2019
License     : GPL-3
Maintainer  : mail@tpapak.com
Stability   : experimental
Portability : POSIX

Following [Rücker's et al. 2013](10.1002/sim.6236) method for substituting a multiarm study with n choose 2 independent reweighted contrasts. The code is the haskell intetrpretation of netmeta's multiarm function
-}

module Data.Meta.Multiarm
  ( reweightMultiArm
  , effectiveVarianceMatrix
  )
where

import           Control.Applicative
import           Data.List
import qualified Data.Map.Strict               as Map
import qualified Data.Set                      as Set
import           Data.Maybe
import           Data.Either
import           Data.Meta.Effects
import qualified Numeric.LinearAlgebra         as LA
import qualified Numeric.LinearAlgebra.Devel   as LAD
import qualified Data.Graph.AdjacencyList      as G
import Data.Meta.Matrices

reweightMultiArm
  :: (Gaussian a, Effect a) => [Contrast a] -> Either String [Contrast a]
reweightMultiArm contrastList =
  if length contrastList > 1 then
    let
      comparisonList = map (\(Contrast ta tb e) -> (ta, tb)) contrastList
      treatments     = foldl' (\vts (ta, tb) -> Set.insert tb (Set.insert ta vts))
                              Set.empty
                              comparisonList
      tsvs     = Map.fromList $ zip (Set.toList treatments) [1 ..]
      vsts     = Map.fromList $ zip [1 ..] (Set.toList treatments)
      edgeList = map
        (\(ta, tb) -> G.fromTuple
          (fromJust (Map.lookup ta tsvs), fromJust (Map.lookup tb tsvs))
        )
        comparisonList
      multigraph   = G.graphFromEdges edgeList
      completeG    = G.completeGraph $ G.numVertices multigraph
      varianceList = map
        (\(Contrast ta tb e) ->
          let
            ed = G.Edge (fromJust (Map.lookup ta tsvs))
                        (fromJust (Map.lookup tb tsvs))
          in  (ed, variance e)
        )
        contrastList
      reverseVarianceList = map (\(ed, v) -> (G.reverseEdge ed, v)) varianceList
      varianceMap         = Map.fromList $ union varianceList reverseVarianceList
      checkEnoughContrasts =
        isJust $ mapM (flip Map.lookup varianceMap) $ G.edges completeG
      checkMoreContrasts = Map.size varianceMap == 2 * length contrastList
    in
      if checkEnoughContrasts && checkMoreContrasts
        then
          let directedComplete = G.removeReverseEdges completeG
              bMatrix          = createBMatrix directedComplete
              variances        = map (\e -> fromJust $ Map.lookup e varianceMap)
                $ G.edges directedComplete
              erMatrix = effectiveVarianceMatrix variances bMatrix
           in case erMatrix of
                Left err -> Left err
                Right rMatrix ->
                  --Lt <- -0.5 * t(B) %*% B %*% R %*% t(B) %*% B / k^2
                  let lt = (-0.5) *
                        (LA.tr bMatrix)
                          <> bMatrix
                          <> rMatrix
                          <> (LA.tr bMatrix)
                          <> bMatrix / (fromIntegral $ G.numVertices directedComplete)^2
                      ll = LA.pinv lt
                      -- W <- diag(diag(L)) - L
                      wv = LA.diag (LA.takeDiag ll) - ll -- ^ Weight matrix
                      vv = 1 / wv -- ^ new variances
                      newContrasts = map (\(Contrast ta tb ef)
                        -> let taindx = fromJust (Map.lookup ta tsvs)
                               tbindx = fromJust (Map.lookup tb tsvs)
                               newVariance = LA.atIndex vv (taindx - 1, tbindx -1)
                               newef = changeVariance ef newVariance
                            in (Contrast ta tb newef)
                       ) contrastList
                   in Right newContrasts
        else Left $ "wrong contrast number" <> show varianceMap
  else Right contrastList

-- | Distribute the edge variances on a symmetrical matrix R given the variances
-- and the edge vertex adjacency matrix B
effectiveVarianceMatrix
  :: [Double] -> LA.Matrix Double -> Either String (LA.Matrix Double)
effectiveVarianceMatrix variances bMatrix =
  let v                 = LA.diagl variances
      nvars             = length variances
      correctDimensions = nvars == fst (LA.size bMatrix)
   in if correctDimensions
        then
          let lp = LA.tr bMatrix <> v <> bMatrix
          in  Right $ LA.diag (LA.takeDiag lp) - lp
        else Left ("wrong dimensions" <> show variances <> show bMatrix)

max_stars_repo_name: tpapak/meta-analysis

max_issues_repo_name: tpapak/meta-analysis

max_forks_repo_name: tpapak/meta-analysis


text:
{-|
Module      : Approx
Description : Implement Approx for Double, Floats and structures
Copyright   : (c) 2020 Kishaloy Neogi
License     : MIT
Maintainer  : Kishaloy Neogi
Email       : nkishaloy@yahoo.com

The library is created to allow for a easy-to-use reasonable way of emulating approx in Haskell. The codes are all in /pure/ Haskell. The idea is to have a natural mathematical feel in writing code, with operators which just works when working with Double and Float and their composite types like lists, vectors etc.

The __Approx__ module defines 2 operators __@=~@__ and __@/~@__, which are for checking /nearly equal to/ and /not nearly equal to/ respectively.

Both the operators __=~__ and __/~__ are put under the class __Approx__.

At least one of the operators have to be defined and the other gets automatically defined.

The library already defines the functions for some of the basic / common types.

For types where __Eq__ is defined like __Char, Bool, Int, Day, Text__ the approx is simply replaced with __==__.

For __Float__ and __Double__, the following formula is used,

@
if max ( |x|, |y| ) < epsilon_Zero
then True
else
  if |x - y| / max ( |x|, |y| ) < epsilon_Eq
  then True
  else False
@

The motivation for defining Approx for classes for which Eq is also defined is to allow for composite types where both Eq and Approx would be present. For example, the following code evaluates to __@True@__, even though the tuple is of type @(Int,Double,Text,[Int],[[Char]],[Double])@.

@
((2,5.35,"happ",[1,2],["ret","we"],[6.78,3.5])
    :: (Int,Double,Text,[Int],[[Char]],[Double]))

    =~ (2,5.35,"happ",[1,2],["ret","we"],[6.78,3.5])
@

For UTCTime, the approx operator checks for equality to the nearest minute. The following expression evaluates to __@True@__.

@
(parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 15:02:15"
    :: Maybe UTCTime)

    =~ parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 15:01:50"
@

The library also provides approx for Complex and common structures like __List, Boxed and Unboxed Vector, Hashmap, Tuples__ and __Maybe__. For all lists, tuples, hashmaps and vectors, the approximation is checked right down to the elements and the order for lists and vectors are important.

For lists, only finite lists are supported. Any use of infinite lists would cause a runtime error.

There are addtional functions __inRange__, __safeInRange__ and __inTol__, which checks for values within Ranges either /explictily/ defined as in __inRange__ and __safeInRange__ or through tolerances as in __inTol__.

You may see the github repository at <https://github.com/n-kishaloy/approx>

-}

{-# LANGUAGE Strict #-}

module Data.Approx
(
-- *How to use this library
-- |Add @approx@ to build-depends and @import Data.Approx@

-- *Documentation
  Approx (..)
, inRange, safeInRange, inTol

) where

import Data.List (foldl')
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as M
import qualified Data.Vector as V
import Data.Time (Day, UTCTime)
import Data.Time.Clock (diffUTCTime)
import qualified Data.HashMap.Strict as Hm
import Data.Hashable

import Data.Text (Text)
import qualified Data.Complex as Cx
import Data.Complex (Complex ( (:+) ) )

eZerD :: Double
eZerD = 1e-8;   {-# INLINE eZerD  #-}

eEqD :: Double
eEqD   = 1e-7;  {-# INLINE eEqD    #-}

eZerF :: Float
eZerF = 1e-6;   {-# INLINE eZerF   #-}

eEqF :: Float
eEqF   = 1e-5;  {-# INLINE eEqF     #-}

infix 4 =~, /~

-- |The class @Approx@ defines 2 operators __@=~@__ and __@/~@__, which are for checking /nearly equal to/ and /not nearly equal to/ respectively.
class Approx a where
  (=~), (/~) :: a -> a -> Bool

  (=~) x y = not (x /~ y)
  {-# INLINE (=~) #-}
  (/~) x y = not (x =~ y)
  {-# INLINE (/~) #-}

  {-# MINIMAL (=~) | (/~) #-}

instance Approx Day where x =~ y = x == y; {-# INLINE (=~) #-}

instance Approx Char where x =~ y = x == y; {-# INLINE (=~) #-}

instance Approx Bool where x =~ y = x == y; {-# INLINE (=~) #-}

instance Approx Text where x =~ y = x == y; {-# INLINE (=~) #-}

instance Approx Int where x =~ y = x == y; {-# INLINE (=~) #-}

instance Approx Integer where x =~ y = x == y; {-# INLINE (=~) #-}

instance Approx UTCTime where
  x =~ y = (round . (/60.0) . realToFrac $ x `diffUTCTime` y) == 0
  {-# INLINE (=~) #-}

instance Approx a => Approx (Cx.Complex a) where
  (a :+ b) =~ (x :+ y) = (a =~ x) && (b =~ y); {-# INLINE (=~) #-}

instance Approx Float where
  x =~ y = (mx < eZerF) || abs (x-y) / mx < eEqF where mx = max (abs x) (abs y)
  {-# INLINE (=~) #-}

instance Approx Double where
  x =~ y = (mx < eZerD) || abs (x-y) / mx < eEqD where mx = max (abs x) (abs y)
  {-# INLINE (=~) #-}

instance Approx a => Approx (Maybe a) where
  Nothing =~ Nothing  =   True
  Just x  =~ Just y   =   x =~ y
  _       =~ _        =   False
  {-# INLINE (=~) #-}

instance Approx a => Approx [a] where
  x =~ y = (length x == length y) && and (zipWith (=~) x y)

instance (Approx a, Approx b) => Approx (a, b) where
  (x,y) =~ (a,b) =
    (   (x =~ a)
    &&  (y =~ b)
    )
  {-# INLINE (=~) #-}

instance (Approx a, Approx b, Approx c) => Approx (a, b, c) where
  (x,y,z) =~ (a,b,c) =
    (   (x =~ a)
    &&  (y =~ b)
    &&  (z =~ c)
    )
  {-# INLINE (=~) #-}

instance (Approx a,Approx b,Approx c,Approx d) => Approx (a,b,c,d) where
  (x,y,z,u) =~ (a,b,c,d) =
    (   (x =~ a)
    &&  (y =~ b)
    &&  (z =~ c)
    &&  (u =~ d)
    )
  {-# INLINE (=~) #-}

instance (Approx a,Approx b,Approx c,Approx d, Approx e) => Approx (a,b,c,d,e) where
  (x,y,z,u,v) =~ (a,b,c,d,e) =
    (   (x =~ a)
    &&  (y =~ b)
    &&  (z =~ c)
    &&  (u =~ d)
    &&  (v =~ e)
    )
  {-# INLINE (=~) #-}

instance (Approx a,Approx b,Approx c,Approx d, Approx e, Approx f) => Approx (a,b,c,d,e,f) where
  (x,y,z,u,v,w) =~ (a,b,c,d,e,f) =
    (   (x =~ a)
    &&  (y =~ b)
    &&  (z =~ c)
    &&  (u =~ d)
    &&  (v =~ e)
    &&  (w =~ f)
    )
  {-# INLINE (=~) #-}

instance (Approx a,Approx b,Approx c,Approx d, Approx e, Approx f, Approx g) => Approx (a,b,c,d,e,f,g) where
  (x,y,z,u,v,w,p) =~ (a,b,c,d,e,f,g) =
    (   (x =~ a)
    &&  (y =~ b)
    &&  (z =~ c)
    &&  (u =~ d)
    &&  (v =~ e)
    &&  (w =~ f)
    &&  (p =~ g)
    )
  {-# INLINE (=~) #-}

instance (Approx a,Approx b,Approx c,Approx d, Approx e, Approx f, Approx g, Approx h) => Approx (a,b,c,d,e,f,g,h) where
  (x,y,z,u,v,w,p,q) =~ (a,b,c,d,e,f,g,h) =
    (   (x =~ a)
    &&  (y =~ b)
    &&  (z =~ c)
    &&  (u =~ d)
    &&  (v =~ e)
    &&  (w =~ f)
    &&  (p =~ g)
    &&  (q =~ h)
    )
  {-# INLINE (=~) #-}

instance (Approx a,Approx b,Approx c,Approx d, Approx e, Approx f, Approx g, Approx h, Approx i) => Approx (a,b,c,d,e,f,g,h,i) where
  (x,y,z,u,v,w,p,q,r) =~ (a,b,c,d,e,f,g,h,i) =
    (   (x =~ a)
    &&  (y =~ b)
    &&  (z =~ c)
    &&  (u =~ d)
    &&  (v =~ e)
    &&  (w =~ f)
    &&  (p =~ g)
    &&  (q =~ h)
    &&  (r =~ i)
    )
  {-# INLINE (=~) #-}

instance (Approx a,Approx b,Approx c,Approx d, Approx e, Approx f, Approx g, Approx h, Approx i, Approx j) => Approx (a,b,c,d,e,f,g,h,i,j) where
  (x,y,z,u,v,w,p,q,r,s) =~ (a,b,c,d,e,f,g,h,i,j) =
    (   (x =~ a)
    &&  (y =~ b)
    &&  (z =~ c)
    &&  (u =~ d)
    &&  (v =~ e)
    &&  (w =~ f)
    &&  (p =~ g)
    &&  (q =~ h)
    &&  (r =~ i)
    &&  (s =~ j)
    )
  {-# INLINE (=~) #-}

instance (Approx a,Approx b,Approx c,Approx d, Approx e, Approx f, Approx g, Approx h, Approx i, Approx j, Approx k) => Approx (a,b,c,d,e,f,g,h,i,j,k) where
  (x,y,z,u,v,w,p,q,r,s,t) =~ (a,b,c,d,e,f,g,h,i,j,k) =
    (   (x =~ a)
    &&  (y =~ b)
    &&  (z =~ c)
    &&  (u =~ d)
    &&  (v =~ e)
    &&  (w =~ f)
    &&  (p =~ g)
    &&  (q =~ h)
    &&  (r =~ i)
    &&  (s =~ j)
    &&  (t =~ k)
    )
  {-# INLINE (=~) #-}

instance (M.Unbox a, Approx a) => Approx (U.Vector a) where
  x =~ y = (U.length x == U.length y) && U.and (U.zipWith (=~) x y)

instance (Approx a) => Approx (V.Vector a) where
  x =~ y = (V.length x == V.length y) && V.and (V.zipWith (=~) x y)

instance (Eq a, Hashable a, Approx b) => Approx (Hm.HashMap a b) where
  x =~ y = fz x y && fz y x where
    fz p q = and $ (\(k,v) -> Hm.lookup k p =~ Just v) <$> Hm.toList q

infix 4 `inRange`, `safeInRange`, `inTol`

{-|@inRange (u,v) x = check if x is inside the range (u,v)@

Note: The function __assumes @u < v@__. This is done to ensure speed of operations. Use safeInRange, otherwise.
-}
inRange :: Ord a => (a,a) -> a -> Bool
inRange (u,v) x = (u <= x) && (x <= v)
{-# INLINE inRange #-}

{-|@safeInRange (u,v) x = check if x is inside the range (u,v)@

Note: The function works even if u>v. However, it has addtional checks and is more expensive. Use only if you are not sure that u < v for your use-case.
-}
safeInRange :: Ord a => (a,a) -> a -> Bool
safeInRange (u,v) = inRange $ if u<v then (u,v) else (v,u)
{-# INLINE safeInRange #-}

{-|@inTol t p x = inRange (p - t, p + t) x@

The Function checks if __@x@__ is close to __@p@__ within a tolerance band __@t@__. Please ensure __@t@__ is /positive/ or there would be /incorrect/ results.
-}
inTol :: (Num a, Ord a) => a -> a -> a -> Bool
inTol t a = inRange (a - t, a + t)
{-# INLINE inTol #-}

max_stars_repo_name: n-kishaloy/approx

max_issues_repo_name: n-kishaloy/approx

max_forks_repo_name: n-kishaloy/approx


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

module FinancialTimeseries.Report.Standard where

import Data.Time (UTCTime)

import Data.Bifunctor (bimap)
import Data.Distributive (Distributive, distribute)

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

import qualified System.Random as R

import qualified Text.Blaze.Html5 as H5

import qualified Graphics.Rendering.Chart.Easy as E

import qualified Statistics.Sample.Histogram as Histo

import qualified FinancialTimeseries.Algorithm.MonteCarlo as AMC
import FinancialTimeseries.Algorithm.Evaluate (Profit, long, evaluate, evaluateFraction)
import FinancialTimeseries.Render.Chart (chartWithEquity)
import qualified FinancialTimeseries.Render.HtmlReader as HtmlReader -- (Config, runHtmlReader)
import FinancialTimeseries.Render.Render (display)
import FinancialTimeseries.Render.Statement (statement, currentTime)
import FinancialTimeseries.Statistics.Statistics (yield, tradeStatistics, stats2cdfChart, stats2list)
import FinancialTimeseries.Test (check_timeseries_prop)
import FinancialTimeseries.Type.Fraction (Fraction(..))
import FinancialTimeseries.Type.Histogram (Histogram(..))
import FinancialTimeseries.Type.Labeled (Labeled(..))
import FinancialTimeseries.Type.MonteCarlo (Broom(..))
import FinancialTimeseries.Type.Strategy (Strategy(..))
import FinancialTimeseries.Type.Timeseries (TimeseriesRaw, first, slice)
import qualified FinancialTimeseries.Type.Timeseries as TS
import FinancialTimeseries.Type.Types (StripPrice, stripPrice, Equity(..), partitionInvested)

import FinancialTimeseries.Util.DistributivePair (distributePair)
import FinancialTimeseries.Util.Pretty (Pretty, pretty)

data Config params gen price a = Config {
  now :: UTCTime
  , reportConfig :: HtmlReader.Config
  , monteCarloConfig :: AMC.Config gen a
  , strategy :: Strategy params price a
  }

report ::
  (Profit price, Distributive price, StripPrice price, Pretty params, R.RandomGen gen, Num a, Fractional a, Real a, E.PlotValue a, Show a, Pretty a, TS.Length (TimeseriesRaw price a)) =>
  Config params gen price a -> TimeseriesRaw price a -> H5.Html
report cfg ts =
  let t = (freeStrategy (strategy cfg)) (parameters (strategy cfg)) ts
      lg = long (partitionInvested (slice t))

      mteCrlo = AMC.mc (monteCarloConfig cfg) lg

      yields =
        let k = fmap (Vec.fromList . map (yield . Vec.map snd))
        in fmap (fmap (bimap k k)) lg

      tradeYields =
        let  j zs = (fmap stats2cdfChart zs, fmap stats2list zs)
             k = j . fmap ((:[]) . Labeled "Trade Yields" . tradeStatistics)
        in fmap (fmap (snd . fmap k)) yields

      tsCharts =
        let convert = Equity . snd . stripPrice
            k = fmap (flip chartWithEquity [t]) . distribute
        in fmap (bimap k k . distributePair) (evaluate (convert (first ts)) lg)

      histogram =
        let histo zs =
              let n = round (fromIntegral (Vec.length zs) / fromIntegral (5 :: Integer) :: Double)
              in Histogram (Histo.histogram n (Vec.map realToFrac zs) :: (Vector Double, Vector Int))
        in fmap (fmap (fmap histo . snd)) yields

      ms = map (\f -> mteCrlo f (evaluateFraction f)) (AMC.fractions (monteCarloConfig cfg))
      ys = fmap (fmap (fmap snd)) (AMC.timeseriesYields ms)
      rdds = fmap (fmap (fmap snd)) (AMC.relativeDrawdowns ms)

      broom = fmap (fmap (fmap (bimap (fmap Broom) (fmap Broom)))) (mteCrlo (Fraction 1.0) evaluate)

      html = HtmlReader.runHtmlReader (reportConfig cfg) $ mconcat $
        currentTime (now cfg)
        : statement ("Timeseries is well formed: " ++ show (check_timeseries_prop t))
        : statement ("Parameters: " ++ pretty (parameters (strategy cfg)))
        : display tsCharts
        : display tradeYields
        : display histogram
        : display broom
        : display ys
        : display rdds
        : []

  in html

max_stars_repo_name: fphh/financial-timeseries

max_issues_repo_name: fphh/financial-timeseries

max_forks_repo_name: fphh/financial-timeseries


text:
module CubicBezier ( CubicBezier (..)
                   , pointAt
                   , bezierSplitAt
                   , isClockwise
                   , inflectionPoints
                   , realInflectionPoint
                   ) where

import Linear
import Control.Lens
import Data.Complex

data CubicBezier = CubicBezier { _p1 :: V2 Double
                               , _c1 :: V2 Double
                               , _c2 :: V2 Double
                               , _p2 :: V2 Double
                               } deriving Show

pointAt :: CubicBezier -> Double -> V2 Double
pointAt bezier t =  ((1 - t) ** 3) *^ _p1 bezier +
                    ((1 - t) ** 2) * 3 * t *^ _c1 bezier +
                    (t ** 2) * (1 - t) * 3 *^ _c2 bezier +
                    (t ** 3) *^ _p2 bezier

bezierSplitAt :: CubicBezier -> Double -> (CubicBezier, CubicBezier)
bezierSplitAt bezier t = (CubicBezier (_p1 bezier) p0 p01 dp, CubicBezier dp p12 p2 (_p2 bezier))
    where
        p0 = _p1 bezier + t *^ (_c1 bezier - _p1 bezier)
        p1 = _c1 bezier + t *^ (_c2 bezier - _c1 bezier)
        p2 = _c2 bezier + t *^ (_p2 bezier - _c2 bezier)

        p01 = p0 + t *^ (p1 - p0)
        p12 = p1 + t *^ (p2 - p1)

        dp = p01 + t *^ (p12 - p01)

isClockwise :: CubicBezier -> Bool
isClockwise bezier = s < 0
    where
        s = (_c1 bezier ^. _x - _p1 bezier  ^. _x) * (_c1 bezier ^. _y + _p1 bezier ^. _y)
          + (_c2 bezier ^. _x - _c1 bezier  ^. _x) * (_c2 bezier ^. _y + _c1 bezier ^. _y)
          + (_p2 bezier ^. _x - _c2 bezier  ^. _x) * (_p2 bezier ^. _y + _c2 bezier ^. _y)
          + (_p1 bezier ^. _x - _p2 bezier  ^. _x) * (_p1 bezier ^. _y + _p2 bezier ^. _y)

inflectionPoints :: CubicBezier -> (Complex Double, Complex Double)
inflectionPoints bezier = (t1, t2)
    where
        pa = _c1 bezier - _p1 bezier
        pb = _c2 bezier - _c1 bezier - pa
        pc = _p2 bezier - _c2 bezier - pa - 2 *^ pb

        a = (pb ^. _x * pc ^. _y - pb ^. _y * pc ^. _x) :+ 0
        b = (pa ^. _x * pc ^. _y - pa ^. _y * pc ^. _x) :+ 0
        c = (pa ^. _x * pb ^. _y - pa ^. _y * pb ^. _x) :+ 0

        t1 = (-b + sqrt (b * b  - 4 * a * c)) / (2 * a)
        t2 = (-b - sqrt (b * b  - 4 * a * c)) / (2 * a)

realInflectionPoint :: Complex Double -> Bool
realInflectionPoint c = imagPart c == 0 && realPart c > 0 && realPart c < 1


text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE StandaloneDeriving    #-}
{-# LANGUAGE TemplateHaskell       #-}
{-# LANGUAGE TypeFamilies          #-}

-- | This module provides type-safe conversions from Foreign C types
--   to Haskell types, and vice versa. It is especially convenient when
--   writing code with Haskell's C FFI.
--
--   This module will not compile on non-64-bit operating systems.
module CCast
  ( Cast(cast)
  , CComplex(CComplex)
  ) where

import Data.Coerce
import Data.Complex
import Data.Int
import Data.Word
import Foreign.C.Types
import Foreign.Ptr
import Foreign.Storable

import Bit64Only

-- | A typeclass for typesafe c-style casts from 'a' to 'b',
--   for use with Haskell's C FFI. This is just a convenience,
--   to prevent from having to think about numeric conversions.
--   Type inference is still weak here; I recommend two solutions:
--
--   With -XTypeApplications:
--
-- @
-- cast \@\CChar \@\Int8 x -- convert x from CChar to Int8
-- @
--
--   With explicit type annotations:
--
-- @
-- cast x :: Int8      -- convert x from CChar to Int8
-- @
--
class Cast a b where
  cast :: Coercible a b => a -> b
  cast = coerce;

-- | Complex number for FFI with the same memory layout as std::complex\<T\>
data CComplex a = CComplex !a !a
  deriving (Eq)

instance Storable a => Storable (CComplex a) where
    sizeOf _ = sizeOf (undefined :: a) * 2
    alignment _ = alignment (undefined :: a)
    poke p (CComplex x y) = do
        pokeElemOff (castPtr p) 0 x
        pokeElemOff (castPtr p) 1 y
    peek p = CComplex
        <$> peekElemOff (castPtr p) 0
        <*> peekElemOff (castPtr p) 1

instance Cast CChar   Int8    where
instance Cast Int8    CChar   where
instance Cast CSChar  Int8    where
instance Cast Int8    CSChar  where
instance Cast CUChar  Word8   where
instance Cast Word8   CUChar  where
instance Cast CShort  Int16   where
instance Cast Int16   CShort  where
instance Cast CUShort Word16  where
instance Cast Word16  CUShort where
instance Cast CInt    Int32   where
instance Cast Int32   CInt    where
instance Cast CUInt   Word32  where
instance Cast Word32  CUInt   where
instance Cast CInt    Int     where; cast = fromIntegral;
-- | WARNING! This conversion is potentially lossy.
instance Cast Int     CInt    where; cast = fromIntegral;
instance Cast CLong   Int64   where
instance Cast Int64   CLong   where
instance Cast CULong  Word64  where
instance Cast Word64  CULong  where
instance Cast CFloat  Float   where
instance Cast Float   CFloat  where
instance Cast CDouble Double  where
instance Cast Double  CDouble where

instance Cast (CComplex CFloat)  (Complex Float)    where; cast (CComplex x y) = cast x :+ cast y; {-# INLINE cast #-}
instance Cast (Complex Float)    (CComplex CFloat)  where; cast (x :+ y) = CComplex (cast x) (cast y); {-# INLINE cast #-}
instance Cast (CComplex CDouble) (Complex Double)   where; cast (CComplex x y) = cast x :+ cast y; {-# INLINE cast #-}
instance Cast (Complex Double)   (CComplex CDouble) where; cast (x :+ y) = CComplex (cast x) (cast y); {-# INLINE cast #-}

instance Cast CPtrdiff   Int64      where
instance Cast Int64      CPtrdiff   where
instance Cast CSize      Word64     where
instance Cast Word64     CSize      where
instance Cast CWchar     Int32      where
instance Cast Int32      CWchar     where
instance Cast CSigAtomic Int32      where
instance Cast Int32      CSigAtomic where
instance Cast CLLong     Int64      where
instance Cast CULLong    Word64     where
instance Cast CBool      Word8      where
instance Cast Word8      CBool      where
instance Cast CIntPtr    Int64      where
instance Cast Int64      CIntPtr    where
instance Cast CUIntPtr   Word64     where
instance Cast Word64     CUIntPtr   where
instance Cast CIntMax    Int64      where
instance Cast Int64      CIntMax    where
instance Cast CClock     Int64      where
instance Cast Int64      CClock     where
instance Cast CTime      Int64      where
instance Cast Int64      CTime      where
instance Cast CUSeconds  Word32     where
instance Cast Word32     CUSeconds  where
instance Cast CSUSeconds Int64      where
instance Cast Int64      CSUSeconds where

-- Proof that we are on a 64-bit architecture.
-- This module will not compile outside of a 64-bit architecture.
myArchitecture :: Refined MustBe64 Int
myArchitecture = $$(refineTH (sizeOf (undefined :: Int)))

max_stars_repo_name: chessai/ccast

max_issues_repo_name: chessai/ccast

max_forks_repo_name: chessai/ccast


text:
-- | Module for providing various functionals used in calculations involving
-- quantum mathematics. This includes the fidelity and trace norm.
module Hoqus.Fidelity where

import Numeric.LinearAlgebra.Data
import Numeric.LinearAlgebra

import Hoqus.MtxFun

-- | Function 'fidelity' calculates the fidelity between two quantum states (or
-- any two square matrices)
--fidelity :: Matrix C -> Matrix C -> C
--fidelity a b = (sum $ map sqrt $ toList $ eigenvalues $ a <> b ) ** 2

-- | Function 'superFidelity' calculates an upper bound for the fidelity using
-- only trace of the products of input matrices.
superFidelity :: Matrix C -> Matrix C -> C
superFidelity a b = (trace (a <> b)) + (sqrt (1 - trace (a <> a))) * (sqrt (1 - trace (b <> b)))

-- | Function 'subFidelity' calculates a lower bound for the fidelity using only
-- race of the products of input matrices.
subFidelity :: Matrix C -> Matrix C -> C
subFidelity a b = (trace p) +  (sqrt 2) * (sqrt ((trace p)*(trace p) - trace (p <> p)))
  where p = a <> b

max_stars_repo_name: jmiszczak/hoqus

max_issues_repo_name: jmiszczak/hoqus

max_forks_repo_name: jmiszczak/hoqus


text:
module Fractal.Control.Controller(
  ApplicationType (UI_Application, Console_Application),
  mandelApplication,
  getApplicationType
)where

import  System.IO (BufferMode (NoBuffering), hSetBuffering, stdout)
import qualified Fractal.Parallel.Evaluator as P
import qualified Fractal.Gui.FractalWindow as W
import qualified Fractal.Utils.Field as F
import Fractal.Math.Mandelbrot (mandel, mandel')
import Data.Complex
import Fractal.Plot.Plotter (imageFractal)
import Data.DateTime

data ApplicationType = UI_Application | Console_Application | None | Empty
  deriving (Show, Eq)

type Point = Complex Double

mandelApplication :: ApplicationType -> IO ()
mandelApplication appType = case appType of
  UI_Application        -> uiApplication
  Console_Application   -> consoleApplication
  Empty                 -> readApplicationTypeFromInput
  None                  -> error "Cant use"

getApplicationType :: [String] -> ApplicationType
getApplicationType [] = Empty
getApplicationType (x:_) = case x of
  "UI_Application"      -> UI_Application
  "Console_Application" -> Console_Application
  _                     -> Empty

readApplicationTypeFromInput :: IO ()
readApplicationTypeFromInput = do
  hSetBuffering stdout NoBuffering
  putStr "Provide application type (UI_Application or Console_Application): "
  applicationType <- getLine
  let appType = getApplicationType [applicationType]
  if appType == Empty then mandelApplication None else mandelApplication appType

consoleApplication :: IO ()
consoleApplication = do
  putStrLn "Running in console mode..."
  input <- consoleReadInput
  coreComputation input
  return ()

coreComputation :: (Double, Int, Double, Point, Point, String) -> IO ()
coreComputation (limit, maxIter, stepSize, sPoint, ePoint, strat) = do
  --let field   = F.generateField sPoint ePoint stepSize
  --let field'  = concat field
  let field'  = F.generateRowField sPoint ePoint stepSize
  let res     = P.mandelEval (P.stringToStrat strat) (mandel' maxIter limit) field'
  let len     = F.getFieldLength sPoint ePoint stepSize
  dateTime    <- getCurrentTime
  let (year, month, day) = toGregorian' dateTime
  imageFractal len res ("./Mandel_" ++ (show year) ++ "_" ++ (show month) ++ "_" ++ (show day) ++ "_" ++ (show $ toSeconds dateTime))
  putStrLn "Image created!"
  return ()

uiApplication :: IO ()
uiApplication = do
  putStrLn "Running in ui mode"
  W.mainW coreComputation
  return ()

consoleReadInput :: IO (Double, Int, Double, Point, Point, String)
consoleReadInput = do
  hSetBuffering stdout NoBuffering
  putStr "Type in recursive limit: "
  fractalLimit <- getLine
  let limit = (read fractalLimit) :: Double
  putStr "Type in max number of iterations: "
  maxIter <- getLine
  let maxIt = (read maxIter) :: Int
  putStr "Type in stepsize: "
  stepSize <- getLine
  let step = (read stepSize) :: Double
  putStr "Type in starting point: "
  startPoint <- getLine
  let startP = (read startPoint) :: Point
  putStr "Type in end point: "
  endPoint <- getLine
  let endP = (read endPoint) :: Point
  putStr "Type in evaluation strategie: "
  evalStrat <- getLine
  return (limit, maxIt, step, startP, endP, evalStrat)

max_stars_repo_name: sigmaticsMUC/FractalHaskell

max_issues_repo_name: sigmaticsMUC/FractalHaskell

max_forks_repo_name: sigmaticsMUC/FractalHaskell


text:
module LearnHaskDA.Chapter5 where
--import Math.Combinatorics.Exact.Binomial
import Numeric.SpecFunctions (choose)
import LearnHaskDA.Chapter2
import LearnHaskDA.Chapter4
import Graphics.EasyPlot
import System.Random
import Data.Number.Erf

{- Using Math.Combinatorics.Exact.Binomial in the exact-combinatorics
package we could have a more general type here and an exact, yet
efficient, computation of choose. But exact-combinatorics is not on
stackage so I changed to the approximation (from Numeric.SpecFunctions
which is on stackage). We only use it for plotting anyway.

-- probabilityMassFunction :: (Integral i, Num n) => i -> i -> n -> n
-}
probabilityMassFunction :: Int -> Int -> Double -> Double
probabilityMassFunction k n p = (n `choose` k)
                              * p^k * (1-p)^(n-k)
plot81 :: IO Bool
plot81 = plot (PNG "coinflips.png") $
              Function2D [Title "Coin Flip Probabilities"]
                         [Range 0 1000]
                         (\k -> probabilityMassFunction (floor k) 1000 0.5)

-- almost equal would be enough
sanitycheck :: Bool
sanitycheck = 1 == sum (map (\k -> probabilityMassFunction k 1000 0.5) [0..1000])

-- sumMid :: (Fractional a, Integral i) => i -> a
sumMid :: Int -> Double
sumMid d = sum $ map (\k -> probabilityMassFunction k 1000 0.5) [500-d .. 500+d]

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

standardDeviation :: (Floating a, Real a) => [a] -> a
standardDeviation values = sqrt (sum [(x-mu)^2 | x <- values]) / sqrt_nm1
  where  mu = average values
         sqrt_nm1 = sqrt $ fromIntegral (length values - 1)

select place = "SELECT "++place++"team, SUM("++place++"score) "++
                 "FROM winloss "++
                 "GROUP BY "++ place ++"team " ++
                 "ORDER BY "++ place ++"team"

test91 = do
  runsAtHome <- queryDatabase "winloss.sql" (select "home")
  runsAway   <- queryDatabase "winloss.sql" (select "away")
  let runsHomeAway = zip (readDoubleColumn runsAtHome 1) (readDoubleColumn runsAway 1)
  plot (PNG "HomeScoreAwayScore.png") $ Data2D [Title "Runs at Home (x axis) and Runs Away (y axis)"] [] runsHomeAway
  let runsHomeAwayDiff = map (\(a,b) -> a-b) runsHomeAway
  plot (PNG "HomeScoreAwayScoreDiff.png") $ Data2D [Title "Difference in Runs at Home and Runs Away"] [] $ zip [1..] runsHomeAwayDiff
  let n = fromIntegral (length runsHomeAwayDiff)
      avg = average runsHomeAwayDiff
      stdDev = standardDeviation runsHomeAwayDiff
      stdErr = stdDev / sqrt n
  return (avg, stdDev, stdErr)

--

from = invnormcdf ((1-95/100)/2 :: Double)
to = - from

test95 = test91 >>= \ (avg, stdDev, stdErr) ->
         return (avg + from*stdErr, avg + to*stdErr)

-- The 95% confidence interval is (-14.9, 26.6) which contains zero so
-- the null hypothesis could not be rejected.

max_stars_repo_name: patrikja/LearningHaskellDataAnalysis

max_issues_repo_name: patrikja/LearningHaskellDataAnalysis

max_forks_repo_name: patrikja/LearningHaskellDataAnalysis


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE LambdaCase #-}

module Year2020.Day20 where

import Util hiding ((??))

import Numeric.LinearAlgebra

import Data.Char
import Data.Function
import Data.List hiding (find)
import Data.List.Extra hiding (find)

import Control.Arrow
import Control.Monad

import System.FilePath

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

type Tile = Matrix Z

makeTile :: [String] -> (Int, Tile)
makeTile (x:xs) = (n, fromLists (map (map (\case '#' -> 1 ; _ -> 0)) xs))
  where n = read (filter isDigit x) :: Int

transforms :: Tile -> [Tile]
transforms m = take 8 . scanl (&) m $ cycle [tr, flipud]

place :: Int -> Map Int Tile -> Map (Int, Int) (Int, Tile)
       -> Int -> Int -> [[[(Int, Tile)]]]
place size tiles placed _ _ | Map.null tiles = pure . transpose
  .  map (map snd) . groupOn (fst . fst) $ Map.toList placed
place size tiles placed x y = do
  let xp = snd (placed Map.! (x - 1, y)) ?? (All, TakeLast 1)
      yp = snd (placed Map.! (x, y - 1)) ?? (TakeLast 1, All)
  (num, tile') <- Map.toList tiles
  tile <- transforms tile'
  guard $ (x == 1 || tile ?? (All, Take 1) == xp)
       && (y == 1 || tile ?? (Take 1, All) == yp)
  let (x', y') = ((+ x) *** (+ 1)) (divMod y size)
  place size (Map.delete num tiles) (Map.insert (x, y) (num, tile) placed) x' y'

part1, part2 :: [[(Int, Tile)]] -> Int
part1 xs = product . map (fst . ($ xs)) $ (.) <$> [head, last] <*> [head, last]
part2 = fromIntegral . minimum . checks . transforms . fromBlocks
  . map (map ((?? (DropLast 1, DropLast 1)) . (?? (Drop 1, Drop 1)) . snd))
  where goal = fromIntegral $ sumElements monster
        check xs = sumElements $ cond (corr2 monster xs) goal 0 goal 0
        checks (x:xs) = map ((sumElements x -) . check) (x:xs)
        monster = snd . makeTile $ "" :
          [ "                  # "
          , "#    ##    ##    ###"
          , " #  #  #  #  #  #   " ]

main = do
  input <- readFile (replaceExtension __FILE__ ".in")
  let tiles = map makeTile . filter notNull $ paragraphs input
      size = round . sqrt . fromIntegral $ length tiles
  print . part2 . head $ place size (Map.fromList tiles) mempty 1 1