text:
--------------------------------------------------------------------------------

{-# OPTIONS_GHC -fno-warn-orphans #-}

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

{-# LANGUAGE AllowAmbiguousTypes       #-}
{-# LANGUAGE BangPatterns              #-}
{-# LANGUAGE CPP                       #-}
{-# LANGUAGE EmptyDataDecls            #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE ForeignFunctionInterface  #-}
{-# LANGUAGE FunctionalDependencies    #-}
{-# LANGUAGE GADTs                     #-}
{-# LANGUAGE KindSignatures            #-}
{-# LANGUAGE LambdaCase                #-}
{-# LANGUAGE MagicHash                 #-}
{-# LANGUAGE MultiParamTypeClasses     #-}
{-# LANGUAGE Rank2Types                #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# LANGUAGE StandaloneDeriving        #-}
{-# LANGUAGE TypeApplications          #-}
{-# LANGUAGE TypeFamilyDependencies    #-}
{-# LANGUAGE TypeInType                #-}
{-# LANGUAGE UnboxedTuples             #-}
{-# LANGUAGE UndecidableInstances      #-}

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

-- | Internal module to Eigen.
--   Here we define all foreign function calls,
--   and some typeclasses integral to the public and private interfaces
--   of the library.
module Eigen.Internal where --   FIXME: Explicit export list

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

import           Control.Monad            (when)
import           Data.Binary              (Binary(put,get))
import           Data.Binary.Get          (getByteString, getWord32be)
import           Data.Binary.Put          (putByteString, putWord32be)
import           Data.Bits                (xor)
import           Data.Complex             (Complex((:+)))
import           Data.Kind                (Type)
import           Data.Proxy               (Proxy(Proxy))
import           Foreign.C.String         (CString, peekCString)
import           Foreign.C.Types          (CInt(CInt), CFloat(CFloat), CDouble(CDouble), CChar)
import           Foreign.ForeignPtr       (ForeignPtr, castForeignPtr, withForeignPtr)
import           Foreign.Ptr              (Ptr, castPtr, nullPtr, plusPtr)
import           Foreign.Storable         (Storable(sizeOf, alignment, poke, peek, peekByteOff, peekElemOff, pokeByteOff, pokeElemOff))
import           GHC.TypeLits             (natVal, KnownNat, Nat)
import           System.IO.Unsafe         (unsafeDupablePerformIO)
import qualified Data.Vector.Storable     as VS
import qualified Data.ByteString          as BS
import qualified Data.ByteString.Internal as BSI

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

-- | Like 'Proxy', but specialised to 'Nat'.
data Row (r :: Nat) = Row
-- | Like 'Proxy', but specialised to 'Nat'.
data Col (c :: Nat) = Col

-- | Used internally. Given a 'KnownNat' constraint, turn the type-level 'Nat' into an 'Int'.
natToInt :: forall n. KnownNat n => Int
{-# INLINE natToInt #-}
natToInt = fromIntegral (natVal @n Proxy)

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

-- | Cast to and from a C-FFI type
--   'Cast' is a closed typeclass with an associated injective type family.
--   It is closed in the sense that we provide only four types
--   with instances for it; and intend for eigen to only be used
--   with those four types. The injectivity of the type family is
--   then useful for avoiding MPTCs. 'Cast' has two functions; 'toC'
--   and 'fromC', where 'toC' goes from a Haskell type to its associated
--   C type for internal use, with the C FFI, and 'fromC' goes from the
--   associated C type to the Haskell type.
class Cast (a :: Type) where
  type family C a = (result :: Type) | result -> a
  toC   :: a -> C a
  fromC :: C a -> a

instance Cast Int where
  type C Int = CInt
  toC = CInt . fromIntegral
  {-# INLINE toC #-}
  fromC (CInt x) = fromIntegral x
  {-# INLINE fromC #-}

instance Cast Float where
  type C Float = CFloat
  toC = CFloat
  {-# INLINE toC #-}
  fromC (CFloat x) = x
  {-# INLINE fromC #-}

instance Cast Double where
  type C Double = CDouble
  toC = CDouble
  {-# INLINE toC #-}
  fromC (CDouble x) = x
  {-# INLINE fromC #-}

instance Cast a => Cast (Complex a) where
  type C (Complex a) = CComplex (C a)
  toC (a :+ b) = CComplex (toC a) (toC b)
  {-# INLINE toC #-}
  fromC (CComplex a b) = (fromC a) :+ (fromC b)
  {-# INLINE fromC #-}

-- | WARNING! 'toC' is lossy for any Int greater than (maxBound :: Int32)!
instance Cast a => Cast (Int, Int, a) where
  type C (Int, Int, a) = CTriplet a
  {-# INLINE toC #-}
  toC (x, y, z) = CTriplet (toC x) (toC y) (toC z)
  {-# INLINE fromC #-}
  fromC (CTriplet x y z) = (fromC x, fromC y, fromC z)

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

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

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

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

-- | FIXME: Doc
data CTriplet a where
  CTriplet :: Cast a => !CInt -> !CInt -> !(C a) -> CTriplet a

deriving instance (Show a, Show (C a)) => Show (CTriplet a)

instance (Storable a, Elem a) => Storable (CTriplet a) where
    sizeOf _ = sizeOf (undefined :: a) + sizeOf (undefined :: CInt) * 2
    alignment _ = alignment (undefined :: CInt)
    poke p (CTriplet row col val) = do
        pokeElemOff (castPtr p) 0 row
        pokeElemOff (castPtr p) 1 col
        pokeByteOff p (sizeOf (undefined :: CInt) * 2) val
    peek p = CTriplet
        <$> peekElemOff (castPtr p) 0
        <*> peekElemOff (castPtr p) 1
        <*> peekByteOff p (sizeOf (undefined :: CInt) * 2)

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

-- | `Elem` is a closed typeclass that encompasses the properties
--   eigen expects its values to possess, and simplifies the external
--   API quite a bit.
class (Num a, Cast a, Storable a, Storable (C a), Code (C a)) => Elem a

instance Elem Float
instance Elem Double
instance Elem (Complex Float)
instance Elem (Complex Double)

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

-- | Encode a C Type as a CInt
--
--   Hack used in FFI wrapper functions when constructing FFI calls
class Code a where; code :: a -> CInt
instance Code CFloat             where; code _ = 0
instance Code CDouble            where; code _ = 1
instance Code (CComplex CFloat)  where; code _ = 2
instance Code (CComplex CDouble) where; code _ = 3

-- | Hack used in constructing FFI calls.
newtype MagicCode = MagicCode CInt deriving Eq

instance Binary MagicCode where
    put (MagicCode _code) = putWord32be $ fromIntegral _code
    get = MagicCode . fromIntegral <$> getWord32be

-- | Hack used in constructing FFI calls.
magicCode :: Code a => a -> MagicCode
magicCode x = MagicCode (code x `xor` 0x45696730)

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

-- | Machine size of a 'CInt'.
intSize :: Int
intSize = sizeOf (undefined :: CInt)

-- | FIXME: Doc
encodeInt :: CInt -> BS.ByteString
encodeInt x = BSI.unsafeCreate (sizeOf x) $ (`poke` x) . castPtr

-- | FIXME: Doc
decodeInt :: BS.ByteString -> CInt
decodeInt (BSI.PS fp fo fs)
    | fs == sizeOf x = x
    | otherwise = error "decodeInt: wrong buffer size"
    where x = performIO $ withForeignPtr fp $ peek . (`plusPtr` fo)

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

-- | 'Binary' instance for 'Data.Vector.Storable.Mutable.Vector'
instance Storable a => Binary (VS.Vector a) where
    put vs = put (BS.length bs) >> putByteString bs where
        (fp,fs) = VS.unsafeToForeignPtr0 vs
        es = sizeOf (VS.head vs)
        bs = BSI.fromForeignPtr (castForeignPtr fp) 0 (fs * es)

    get = get >>= getByteString >>= \bs -> let
        (fp,fo,fs) = BSI.toForeignPtr bs
        es = sizeOf (VS.head vs)
        -- `plusForeignPtr` is used qualified here to just remind a reader
        -- that it is defined internally within eigen
        vs = VS.unsafeFromForeignPtr0 (Eigen.Internal.plusForeignPtr fp fo) (fs `div` es)
        in return vs

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

-- | FIXME: Doc
data CSparseMatrix a
-- | FIXME: Doc
type CSparseMatrixPtr a = Ptr (CSparseMatrix a)

-- | FIXME: Doc
data CSolver a
-- | FIXME: Doc
type CSolverPtr a = Ptr (CSolver a)

-- {-# INLINE unholyPerformIO #-}
-- unholyPerformIO :: IO a -> a
-- unholyPerformIO (IO m) = case m realWorld# of (# _, r #) -> r

-- | FIXME: replace with unholyPerformIO (?)
performIO :: IO a -> a
performIO = unsafeDupablePerformIO

-- | FIXME: Doc
plusForeignPtr :: ForeignPtr a -> Int -> ForeignPtr b
plusForeignPtr fp fo = castForeignPtr fp1 where
    vs :: VS.Vector CChar
    vs = VS.unsafeFromForeignPtr (castForeignPtr fp) fo 0
    (fp1, _) = VS.unsafeToForeignPtr0 vs

foreign import ccall "eigen-proxy.h free" c_freeString :: CString -> IO ()

call :: IO CString -> IO ()
call func = func >>= \c_str -> when (c_str /= nullPtr) $
    peekCString c_str >>= \str -> c_freeString c_str >> fail str

foreign import ccall "eigen-proxy.h free" free :: Ptr a -> IO ()

foreign import ccall "eigen-proxy.h eigen_setNbThreads" c_setNbThreads :: CInt -> IO ()
foreign import ccall "eigen-proxy.h eigen_getNbThreads" c_getNbThreads :: IO CInt

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

#let api1 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall a . Code (C a) => %s\n%s = c_%s (code (undefined :: (C a)))", #name, #name, args, #name, args, #name, #name

#api1 random,        "Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 identity,      "Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 add,           "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 sub,           "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 mul,           "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 diagonal,      "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 transpose,     "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 inverse,       "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 adjoint,       "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 conjugate,     "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 normalize,     "Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 sum,           "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 prod,          "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 mean,          "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 norm,          "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 trace,         "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 squaredNorm,   "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 blueNorm,      "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 hypotNorm,     "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 determinant,   "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 rank,          "CInt -> Ptr CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 image,         "CInt -> Ptr (Ptr (C a)) -> Ptr CInt -> Ptr CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 kernel,        "CInt -> Ptr (Ptr (C a)) -> Ptr CInt -> Ptr CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 solve,         "CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api1 relativeError, "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"

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

#let api2 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall a . Code (C a) => %s\n%s = c_%s (code (undefined :: (C a)))", #name, #name, args, #name, args, #name, #name

#api2 sparse_new,           "CInt -> CInt -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_clone,         "CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_fromList,      "CInt -> CInt -> Ptr (CTriplet a) -> CInt -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_toList,        "CSparseMatrixPtr a -> Ptr (CTriplet a) -> CInt -> IO CString"
#api2 sparse_free,          "CSparseMatrixPtr a -> IO CString"
#api2 sparse_makeCompressed,"CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_uncompress,    "CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_isCompressed,  "CSparseMatrixPtr a -> Ptr CInt -> IO CString"
#api2 sparse_transpose,     "CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_adjoint,       "CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_pruned,        "CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_prunedRef,     "CSparseMatrixPtr a -> Ptr (C a) -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_scale,         "CSparseMatrixPtr a -> Ptr (C a) -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_nonZeros,      "CSparseMatrixPtr a -> Ptr CInt -> IO CString"
#api2 sparse_innerSize,     "CSparseMatrixPtr a -> Ptr CInt -> IO CString"
#api2 sparse_outerSize,     "CSparseMatrixPtr a -> Ptr CInt -> IO CString"
#api2 sparse_coeff,         "CSparseMatrixPtr a -> CInt -> CInt -> Ptr (C a) -> IO CString"
#api2 sparse_coeffRef,      "CSparseMatrixPtr a -> CInt -> CInt -> Ptr (Ptr (C a)) -> IO CString"
#api2 sparse_cols,          "CSparseMatrixPtr a -> Ptr CInt -> IO CString"
#api2 sparse_rows,          "CSparseMatrixPtr a -> Ptr CInt -> IO CString"
#api2 sparse_norm,          "CSparseMatrixPtr a -> Ptr (C a) -> IO CString"
#api2 sparse_squaredNorm,   "CSparseMatrixPtr a -> Ptr (C a) -> IO CString"
#api2 sparse_blueNorm,      "CSparseMatrixPtr a -> Ptr (C a) -> IO CString"
#api2 sparse_add,           "CSparseMatrixPtr a -> CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_sub,           "CSparseMatrixPtr a -> CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_mul,           "CSparseMatrixPtr a -> CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_block,         "CSparseMatrixPtr a -> CInt -> CInt -> CInt -> CInt -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_fromMatrix,    "Ptr (C a) -> CInt -> CInt -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api2 sparse_toMatrix,      "CSparseMatrixPtr a -> Ptr (C a) -> CInt -> CInt -> IO CString"
#api2 sparse_values,        "CSparseMatrixPtr a -> Ptr CInt -> Ptr (Ptr (C a)) -> IO CString"
#api2 sparse_outerStarts,   "CSparseMatrixPtr a -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"
#api2 sparse_innerIndices,  "CSparseMatrixPtr a -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"
#api2 sparse_innerNNZs,     "CSparseMatrixPtr a -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"
#api2 sparse_setZero,       "CSparseMatrixPtr a -> IO CString"
#api2 sparse_setIdentity,   "CSparseMatrixPtr a -> IO CString"
#api2 sparse_reserve,       "CSparseMatrixPtr a -> CInt -> IO CString"
#api2 sparse_resize,        "CSparseMatrixPtr a -> CInt -> CInt -> IO CString"

#api2 sparse_conservativeResize,    "CSparseMatrixPtr a -> CInt -> CInt -> IO CString"
#api2 sparse_compressInplace,       "CSparseMatrixPtr a -> IO CString"
#api2 sparse_uncompressInplace,     "CSparseMatrixPtr a -> IO CString"

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

#let api3 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> CInt -> %s\n%s :: forall s a . (Code s, Code (C a)) => s -> %s\n%s s = c_%s (code (undefined :: (C a))) (code s)", #name, #name, args, #name, args, #name, #name

#api3 sparse_la_newSolver,          "Ptr (CSolverPtr a) -> IO CString"
#api3 sparse_la_freeSolver,         "CSolverPtr a -> IO CString"
#api3 sparse_la_factorize,          "CSolverPtr a -> CSparseMatrixPtr a -> IO CString"
#api3 sparse_la_analyzePattern,     "CSolverPtr a -> CSparseMatrixPtr a -> IO CString"
#api3 sparse_la_compute,            "CSolverPtr a -> CSparseMatrixPtr a -> IO CString"
#api3 sparse_la_tolerance,          "CSolverPtr a -> Ptr CDouble -> IO CString"
#api3 sparse_la_setTolerance,       "CSolverPtr a -> CDouble -> IO CString"
#api3 sparse_la_maxIterations,      "CSolverPtr a -> Ptr CInt -> IO CString"
#api3 sparse_la_setMaxIterations,   "CSolverPtr a -> CInt -> IO CString"
#api3 sparse_la_info,               "CSolverPtr a -> Ptr CInt -> IO CString"
#api3 sparse_la_error,              "CSolverPtr a -> Ptr CDouble -> IO CString"
#api3 sparse_la_iterations,         "CSolverPtr a -> Ptr CInt -> IO CString"
#api3 sparse_la_solve,              "CSolverPtr a -> CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
-- #api3 sparse_la_solveWithGuess,     "CSolverPtr a -> CSparseMatrixPtr a -> CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api3 sparse_la_matrixQ,            "CSolverPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api3 sparse_la_matrixR,            "CSolverPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api3 sparse_la_setPivotThreshold,  "CSolverPtr a -> CDouble -> IO CString"
#api3 sparse_la_rank,               "CSolverPtr a -> Ptr CInt -> IO CString"
#api3 sparse_la_matrixL,            "CSolverPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api3 sparse_la_matrixU,            "CSolverPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"
#api3 sparse_la_setSymmetric,       "CSolverPtr a -> CInt -> IO CString"
#api3 sparse_la_determinant,        "CSolverPtr a -> Ptr (C a) -> IO CString"
#api3 sparse_la_logAbsDeterminant,  "CSolverPtr a -> Ptr (C a) -> IO CString"
#api3 sparse_la_absDeterminant,     "CSolverPtr a -> Ptr (C a) -> IO CString"
#api3 sparse_la_signDeterminant,    "CSolverPtr a -> Ptr (C a) -> IO CString"

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

ext: hsc

max_stars_repo_name: nilsalex/eigen

max_issues_repo_name: nilsalex/eigen

max_forks_repo_name: nilsalex/eigen


text:
module Object.Vertex (Vertex (..), createVertex) where

import Numeric.Vector
import Types

createVertex :: Vec3d -> Vec3d -> Vec2d -> Vec3d -> Vec3d -> Vertex
createVertex pos color texcoord normal tangent = Vertex { v_getPos = pos
                                                        , v_getColor = color
                                                        , v_getTexCoord = texcoord
                                                        , v_getNormal = normal
                                                        , v_getTangent = tangent
                                                        , v_getBitangent = tangent × normal
                                                        }

ext: hs


text:
module FinancialTimeseries.Algorithm.MovingAverage where

import qualified Data.Vector as Vec

import qualified Statistics.Sample as Sample

import FinancialTimeseries.Type.Labeled (Labeled(..))
import FinancialTimeseries.Type.Segment (Segment(..), HalfSegment(..))
import FinancialTimeseries.Type.Strategy (Strategy(..))
import FinancialTimeseries.Type.Timeseries (TimeseriesRaw(..), Timeseries(..))
import FinancialTimeseries.Type.Types (StripPrice, stripPrice)

import FinancialTimeseries.Util.Pretty (Pretty, pretty)
import FinancialTimeseries.Util.ToFileString (ToFileString, toFileString)

newtype Window = Window {
  unWindow :: Int
  } deriving (Show)

instance ToFileString Window where
  toFileString (Window w) = "Window-" ++ show w

instance Pretty Window where
  pretty = show

data UpDown =
  Up Int
  | Down Int
  deriving (Show)

movingAverage ::
  (StripPrice price, Real a, Fractional a) =>
  Window -> TimeseriesRaw price a -> Timeseries price a
movingAverage (Window m) ts@(TimeseriesRaw _ pvs) =
  let vs = stripPrice pvs
      idx = Vec.fromList [0 .. Vec.length vs - m]

      slc i =
        let ws = Vec.slice i m vs
            (tn, xn) = Vec.last ws
        in (tn, realToFrac (Sample.mean (Vec.map (realToFrac . snd) ws)), xn)
      us = Vec.map slc idx

      f i ((_, m0, x0), (_, m1, x1)) acc
        | m0 <= x0 && m1 > x1 = Up (i+m) : acc
        | m0 > x0 && m1 <= x1 = Down (i+m) : acc
        | otherwise = acc

      xs = Vec.ifoldr' f [] (Vec.zip us (Vec.tail us))

      g [] = (Nothing, [])
      g [Down i] = (Just (HalfSegment i), [])
      g (Down i:Up j:zs) = fmap (Segment i j :) (g zs)
      g (Up _:zs) = g zs
      g (Down _:zs) = g zs

      h (t, s, _) = (t, s)

      (hs, isegs) = g xs

      res = Timeseries {
        timeseriesRaw = ts
        , investedSegments = isegs
        , lastSegment = hs
        , additionalSeries = [Labeled "Moving Average" (Vec.map h us)]
        }
  in res

movingAverageStrategy ::
  (StripPrice price, Real a, Fractional a) =>
  Window -> Strategy Window price a
movingAverageStrategy w = Strategy "MovingAverage" w movingAverage

max_stars_repo_name: fphh/financial-timeseries

max_issues_repo_name: fphh/financial-timeseries

max_forks_repo_name: fphh/financial-timeseries


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE OverloadedLabels #-}

module NotificationCenter where

import NotificationCenter.Notification
  (DisplayingNotificaton(..), showNotification, updateNoti)
import NotificationCenter.Notifications
  (NotifyState(..), startNotificationDaemon, Notification(..)
  , hideAllNotis)
import NotificationCenter.Glade (glade, style)
import NotificationCenter.Button
  (Button(..), createButton, setButtonState)
import TransparentWindow
import Helpers
import NotificationCenter.Notifications.Data

import Prelude

import Data.Int (Int32(..))
import Data.Tuple.Sequence (sequenceT)
import Data.Maybe
import Data.IORef
import Data.List
import Data.List.Split (splitOn)
import Data.Time
import Data.Time.LocalTime
import qualified Data.Text as Text
import qualified Data.ByteString.Char8 as BS
import qualified Data.Map as Map
import Data.Complex
import Data.Monoid ((<>))

import Control.Monad
import Control.Applicative
import Control.Concurrent (forkIO, threadDelay, ThreadId(..))
import Control.Concurrent.STM
  ( readTVarIO, modifyTVar', TVar(..), atomically, newTVarIO )

import System.Process (runCommand)
import System.Locale.Read
import System.Posix.Signals (sigUSR1)
import System.Posix.Daemonize (serviced, daemonize)
import System.Directory (getHomeDirectory)

import DBus ( fromVariant )

import GI.Gtk
       (widgetSetHalign, widgetSetHexpand, buttonNew, setWidgetMargin, buttonSetRelief, widgetSetSizeRequest, widgetShowAll, widgetShow, widgetHide, onWidgetDestroy
       , windowSetDefaultSize, setWindowTitle, boxPackStart, boxNew
       , setWindowWindowPosition, WindowPosition(..), windowMove
       , frameSetShadowType, aspectFrameNew
       , widgetGetAllocatedHeight, widgetGetAllocatedWidth, onWidgetDraw
       , onWidgetLeaveNotifyEvent, onWidgetMotionNotifyEvent
       , widgetAddEvents, alignmentSetPadding, alignmentNew, rangeSetValue
       , scaleSetDigits, scaleSetValuePos, rangeGetValue
       , afterScaleButtonValueChanged, scaleNewWithRange, containerAdd
       , buttonBoxNew, mainQuit, onButtonActivate
       , toggleButtonGetActive, onToggleButtonToggled, buttonSetUseStock
       , toggleButtonNewWithLabel, onButtonClicked
       , buttonNewWithLabel, widgetQueueDraw, drawingAreaNew
       , windowNew, widgetDestroy, dialogRun, setAboutDialogComments
       , setAboutDialogAuthors, setAboutDialogVersion
       , setAboutDialogProgramName, aboutDialogNew, labelNew, get
       , afterWindowSetFocus, labelSetText
       , onWidgetFocusOutEvent, onWidgetKeyReleaseEvent, widgetGetParentWindow
       , onButtonClicked, windowGetScreen, boxNew, widgetSetValign)
import qualified GI.Gtk as Gtk (containerAdd, Window(..), Box(..), Label(..), Button(..))

import qualified GI.Gtk as GI (init, main)
import GI.GLib (sourceRemove, timeoutAdd, unixSignalAdd)
import GI.GLib.Constants
import GI.Gdk.Constants
import GI.Gdk.Flags (EventMask(..))
import GI.Gtk.Enums
       (Orientation(..), WindowType(..), ShadowType(..)
       , PositionType(..), ReliefStyle(..), Align(..))
import Data.GI.Base.BasicConversions (gflagsToWord)
import qualified GI.Gdk.Objects.Window

data State = State
  { stMainWindow :: Gtk.Window
  , stNotiBox :: Gtk.Box
  , stTimeLabel :: Gtk.Label
  , stDateLabel :: Gtk.Label
  , stDeleteAll :: Gtk.Button
  , stUserButtons :: [ Button ]
  , stNotiState :: TVar NotifyState
  , stDisplayingNotiList :: [ DisplayingNotificaton ]
  , stNotisForMe :: [ Notification ]
  , stCenterShown :: Bool
  }

setTime :: TVar State -> IO Bool
setTime tState = do
  state <- readTVarIO tState
  now <- zonedTimeToLocalTime <$> getZonedTime
  zone <- System.Locale.Read.getCurrentLocale
  let format = Text.pack . flip (formatTime zone) now
  labelSetText (stTimeLabel state) $ format "%H:%M"
  labelSetText (stDateLabel state) $ format "%A, %x"
  return False

startSetTimeThread :: TVar State -> IO ()
startSetTimeThread tState = do
  runAfterDelay 1000 (startSetTimeThread' tState)
  return ()

startSetTimeThread' :: TVar State -> IO ()
startSetTimeThread' tState = do
  addSource (setTime tState)
  time <- fromIntegral <$> diffTimeToPicoseconds <$> utctDayTime <$> getCurrentTime
  let delay = (60 * 1000000) - ((ceiling (time / 1000000)) `mod` (1000000 * 60))
  runAfterDelay delay (startSetTimeThread' tState)
  return ()

createNotiCenter :: TVar State -> Config -> IO ()
createNotiCenter tState config = do
  objs <- createTransparentWindow (Text.pack glade)
    [ "main_window"
    , "label_time"
    , "label_date"
    , "box_notis"
    , "box_buttons"
    , "button_deleteAll" ]
    (Just "Notification area")

  mainWindow <- window objs "main_window"
  notiBox <- box objs "box_notis"
  buttonBox <- box objs "box_buttons"
  timeL <- label objs "label_time"
  timeD <- label objs "label_date"
  deleteButton <- button objs "button_deleteAll"

  screen <- windowGetScreen mainWindow
  setStyle screen $ BS.pack $
    replaceColors config style
--  (Just mainWindowGDK) <- widgetGetParentWindow mainWindow

  onButtonClicked deleteButton $ do
    displayList <- stDisplayingNotiList <$> readTVarIO tState
    mapM (removeNoti tState) displayList
    return ()

  let buttons = zip
        (split $ removeOuterLetters $ configLabels config)
        (split $ removeOuterLetters $ configCommands config)
      (buttonLabels, buttonCommands) = unzip buttons
      margin = fromIntegral $ configButtonMargin config
      width = fromIntegral (((configWidth config) - 20)
                            `div` (configButtonsPerRow config))
              - margin * 2
      height = fromIntegral $ configButtonHeight config
      linesNeeded = fromIntegral $ ceiling
        $ ((fromIntegral $ length buttons) / (fromIntegral $ configButtonsPerRow config))
  lines' <- sequence $ take linesNeeded $ repeat $ boxNew OrientationHorizontal 0
  buttons' <- sequence $ map
    (\(label, command) -> createButton config width height command label)
    buttons

  sequence $ map (\(box, buttons'') ->
                    sequence $ Gtk.containerAdd box <$> buttons'')
    $ zip (reverse lines') (splitEvery (configButtonsPerRow config)
                            $ (buttonButton <$> buttons'))
  sequence $ Gtk.containerAdd buttonBox <$> lines'
  sequence $ widgetShowAll <$> lines'

  atomically $ modifyTVar' tState $
    \state -> state { stMainWindow = mainWindow
                      , stNotiBox = notiBox
                      , stTimeLabel = timeL
                      , stDateLabel = timeD
                      , stDeleteAll = deleteButton
                      , stNotisForMe = []
                      , stUserButtons = buttons' }

  startSetTimeThread tState

  onWidgetLeaveNotifyEvent mainWindow $ \(_) -> do
    hideNotiCenter tState
    return True

  --  (screenH, screenW) <- getScreenProportions mainWindow
  (screenW, screenY, screenH) <- getScreenPos mainWindow
    (fromIntegral $ configNotiCenterMonitor config)

  windowSetDefaultSize mainWindow
    (fromIntegral $ configWidth config) -- w
    (screenH - barHeight) -- h
  windowMove mainWindow
    (screenW - (fromIntegral $ configWidth config)) -- x
    (screenY + barHeight)  -- y
  onWidgetDestroy mainWindow mainQuit
  return ()
    where barHeight = fromIntegral $ configBarHeight config

parseNotisForMe tState = do
  state <- readTVarIO tState
  -- do stuff
  let buttons = stUserButtons state
      notisForMe = stNotisForMe state
      myNotiHints = notiHints <$> notisForMe
      maybeIds = Map.lookup "id" <$> myNotiHints
      ids = join <$> map (\id' -> fromVariant <$> id') maybeIds :: [ Maybe Int32 ]
      ids' = map checkId ids
        where
          checkId (Just id)
            | id < (fromIntegral $ length buttons) = Just id
            | otherwise = Nothing
          checkId Nothing = Nothing
      maybeStates = Map.lookup "state" <$> myNotiHints
      states = join <$> map (\id' -> fromVariant <$> id') maybeStates :: [ Maybe Bool ]
      comb = sequenceT <$> (zip
                            (map
                             (\id -> (!!) buttons <$> (fromIntegral <$> id))
                             ids') states) :: [ Maybe (Button, Bool) ]
  sequence $ map (\comb' -> sequence $
                    (\(button, state) ->
                           setButtonState button state) <$> comb') comb
  atomically $ modifyTVar' tState
    (\state -> state { stNotisForMe = [] })
  return False

hideNotiCenter tState = do
  state <- readTVarIO tState
  mainWindow <- stMainWindow <$> readTVarIO tState
  widgetHide mainWindow
  atomically $ modifyTVar' tState
    (\state -> state { stCenterShown = False })

showNotiCenter tState notiState = do
  state <- readTVarIO tState
  mainWindow <- stMainWindow <$> readTVarIO tState
  newShown <- if stCenterShown state then
    do
      widgetHide mainWindow
      return False
    else
    do
      hideAllNotis $ stNotiState state
      widgetShow mainWindow
      return True
  atomically $ modifyTVar' tState
    (\state -> state {stCenterShown = newShown })
  return True

updateNotisForMe :: TVar State -> IO()
updateNotisForMe tState = do
  state <- readTVarIO tState
  notiState <- readTVarIO $ stNotiState state
  atomically $ modifyTVar' (stNotiState state) (
    \state -> state { notiForMeList = [] })
  atomically $ modifyTVar' tState (
    \state -> state { stNotisForMe = notiForMeList notiState
                                     ++ stNotisForMe state})

  addSource (parseNotisForMe tState)
  return ()

updateNotis :: Config -> TVar State -> IO()
updateNotis config tState = do
  state <- readTVarIO tState
  notiState <- readTVarIO $ stNotiState state

  let newNotis = filter (
        \n -> (find (\nd -> dNotiId nd == notiId n )
               (stDisplayingNotiList state))
              == Nothing
              && ((configIgnoreTransient config) || (not $ notiTransient n)))
                 $ notiStList notiState
  newNotis' <- mapM (
    \n -> do
      let newNoti = DisplayingNotificaton
                    { dNotiId = notiId n }
      showNotification config
        (stNotiBox state) newNoti
        (stNotiState state) $ removeNoti tState
    ) newNotis

  let delNotis = filter (\nd -> (find (\n -> dNotiId nd == notiId n)
                                $ notiStList notiState) == Nothing)
                 $ stDisplayingNotiList state
  atomically $ modifyTVar' tState (
    \state -> state { stDisplayingNotiList =
                      newNotis' ++ stDisplayingNotiList state
                    })
  setDeleteAllState tState
  mapM (removeNoti tState) $ delNotis
  when (stCenterShown state) $
    do
      hideAllNotis $ stNotiState state
  return ()

removeNoti :: TVar State -> DisplayingNotificaton -> IO ()
removeNoti tState dNoti = do
  state <- readTVarIO tState
  atomically $ modifyTVar' tState $ \state ->
    state { stDisplayingNotiList =
            filter ((/=) dNoti) $ stDisplayingNotiList state}
  setDeleteAllState tState
  atomically $ modifyTVar' (stNotiState state) $ \state ->
    state { notiStList = filter
                         (\n -> notiId n /= dNotiId dNoti) $
                         notiStList state }
  addSource $ do
    dNotiDestroy dNoti
    return False
  return ()

setDeleteAllState tState = do
  addSource $ do
    state <- readTVarIO tState
    if (length $ stDisplayingNotiList state) > 1 then
      widgetShow $ stDeleteAll state
      else
      widgetHide $ stDeleteAll state
    return False
  return ()

getConfig p =
  Config
  {
    -- notification-center
    configBarHeight = r 0 p nCenter "marginTop"
  , configWidth = r 500 p nCenter "width"
  , configStartupCommand = r' "" p nCenter "startupCommand"
  , configNotiCenterMonitor = r 0 p nCenter "monitor"
  , configNotiCenterNewFirst = r'' True p nCenter "newFirst"
  , configIgnoreTransient = r'' False p nCenter "ignoreTransient"
  , configMatchingRules = zip3 match (modify ++ repeat id) $ run ++ repeat Nothing -- run
  , configNotiMarkup = r'' True p nCenter "useMarkup"
  , configSendNotiClosedDbusMessage = r'' False p nCenter "configSendNotiClosedDbusMessage"

    -- notification-center-notification-popup
  , configNotiDefaultTimeout = r 10000 p nPopup "notiDefaultTimeout"
  , configDistanceTop = r 50 p nPopup "distanceTop"
  , configDistanceRight = r 50 p nPopup "distanceRight"
  , configDistanceBetween = r 20 p nPopup "distanceBetween"
  , configWidthNoti = r 300 p nPopup "width"
  , configNotiMonitor = r 0 p nPopup "monitor"

    -- buttons
  , configButtonsPerRow = r 5 p buttons "buttonsPerRow"
  , configButtonHeight = r 60 p buttons "buttonHeight"
  , configButtonMargin = r 2 p buttons "buttonMargin"
  , configLabels = r' "" p buttons "labels"
  , configCommands = r' "" p buttons "commands"
  , configUserButtonColor = r' "#fee" p buttons "buttonColor"
  , configUserButtonHover = r' "rgba(0, 20, 20, 0.2)" p buttons "buttonHover"
  , configUserButtonHoverColor = r' "#fee" p buttons "buttonHoverColor"
  , configUserButtonBackground = r' "rgba(255, 255, 255, 0.15)" p buttons "buttonBackground"
  , configUserButtonTextSize = r' "12px" p buttons "buttonTextSize"
  , configUserButtonState1 = r' "rgba(255,255,255,0.5)" p buttons "buttonState1"
  , configUserButtonState2 = r' "rgba(255,255,255,0.3)" p buttons "buttonState2"
  , configUserButtonState1Color = r' "#fff" p buttons "buttonState1Color"
  , configUserButtonState2Color = r' "#fff" p buttons "buttonState2Color"
  , configUserButtonState1Hover = r' "rgba(0, 20, 20, 0.4)" p buttons "buttonState1Hover"
  , configUserButtonState2Hover = r' "rgba(0, 20, 20, 0.3)" p buttons "buttonState2Hover"
  , configUserButtonState1HoverColor = r' "#fee" p buttons "buttonState1HoverColor"
  , configUserButtonState2HoverColor = r' "#fee" p buttons "buttonState2HoverColor"

    -- colors
  , configBackground = r' "rgba(29, 27, 20, 0.6)" p colors "background"
  , configBackgroundNoti = r' "rgba(9, 0, 0, 0.5)" p colors "notiBackground"
  , configNotiLabelColor = r' "#fef3f6" p colors "notiColor"
  , configCritical = r' "rgba(255, 0, 0, 0.5)" p colors "critical"
  , configCriticalInCenter = r' "rgba(155, 0, 20, 0.5)" p colors "criticalInCenter"
  , configButtonColor = r' "#FFFFFF" p colors "buttonColor"
  , configButtonHover = r' "rgba(0, 20, 20, 0.2)" p colors "buttonHover"
  , configButtonHoverColor = r' "#fee" p colors "buttonHoverColor"
  , configButtonBackground = r' "transparent" p colors "buttonBackground"
  , configLabelColor = r' "#eae2e0" p colors "labelColor"
  , configCriticalColor = r' "#FFFFFF" p colors "criticalColor"
  , configCriticalInCenterColor = r' "#FFFFFF" p colors "criticalInCenterColor"
  }
  where nPopup = "notification-center-notification-popup"
        nCenter = "notification-center"
        colors = "colors"
        buttons = "buttons"
        r = readConfig
        r' = readConfig
        r'' = readConfig
        keys conditions = map (\c -> if length (splitted c) == 2
                                     then (splitted c) else ["", ""]) conditions
          where splitted condition = splitOn "=" condition

        match = map (matcherFunction) $ splitOn ";" <$>
                (split $ removeOuterLetters $ r' "" p nCenter "match")
                where lookupFun name noti = Text.unpack $ fromMaybe (Text.pack "")
                        ((lookup name
                          [ ("title", notiSummary)
                          , ("body", notiBody)
                          , ("app", notiAppName)
                          , ("time", notiTime)
                          , ("icon", notiIcon) ]) <*> (Just noti))
                      matcherFunction conditions = \noti -> foldl (
                        \matches (k:v:[]) -> matches && ((v == lookupFun k noti)))
                                                      True (keys conditions)
        modify = map (matcherFunction) $ splitOn ";" <$>
                 (split $ removeOuterLetters $ r' "" p nCenter "modify")
          where matcherFunction condition = \noti -> foldl (
                  \noti (k:v:[]) -> switch k v noti) noti (keys condition)
                        where switch k v noti
                                | k == "title" = noti { notiSummary = Text.pack v }
                                | k == "body" = noti { notiBody = Text.pack v }
                                | k == "app" = noti { notiAppName = Text.pack v }
                                | k == "time" = noti { notiTime = Text.pack v }
                                | k == "icon" = noti { notiIcon = Text.pack v }
                                | k == "transient" && v == "true" = noti { notiTransient = True }
                                | k == "transient" && v == "false" = noti { notiTransient = False }
                                | k == "noClosedMsg" && v == "true" = noti { notiSendClosedMsg = False }
                                | otherwise = noti

        run = [Nothing]

replaceColors config style =
  replace "replaceme0000" (configBackground config) $
  replace "replaceme0001" (configCritical config) $
  replace "replaceme0002" (configCriticalInCenter config) $
  replace "replaceme0003" (configButtonColor config) $
  replace "replaceme0004" (configButtonHover config) $
  replace "replaceme0005" (configButtonHoverColor config) $
  replace "replaceme0006" (configButtonBackground config) $
  replace "replaceme0007" (configLabelColor config) $
  replace "replaceme0008" (configCriticalColor config) $
  replace "replaceme0009" (configCriticalInCenterColor config) $
  replace "replaceme0010" (configBackgroundNoti config) $
  replace "replaceme0011" (configNotiLabelColor config) $
  replace "replaceme0012" (configUserButtonTextSize config) $
  replace "replaceme0013" (configUserButtonBackground config) $
  replace "replaceme0014" (configUserButtonColor config) $
  replace "replaceme0015" (configUserButtonHover config) $
  replace "replaceme0016" (configUserButtonHoverColor config) $
  replace "replaceme0017" (configUserButtonState1 config) $
  replace "replaceme0018" (configUserButtonState2 config) $
  replace "replaceme0019" (configUserButtonState1Color config) $
  replace "replaceme0020" (configUserButtonState2Color config) $
  replace "replaceme0021" (configUserButtonState1Hover config) $
  replace "replaceme0022" (configUserButtonState2Hover config) $
  replace "replaceme0023" (configUserButtonState1HoverColor config) $
  replace "replaceme0024" (configUserButtonState2HoverColor config) style

getInitialState = do
  newTVarIO $ State
    { stDisplayingNotiList = []
    , stCenterShown = False}

main' :: IO ()
main' = do
  GI.init Nothing

  homeDir <- getHomeDirectory
  config <- getConfig <$> (readConfigFile
    (homeDir ++ "/.config/deadd/deadd.conf"))

  istate <- getInitialState
  notiState <- startNotificationDaemon config
    (updateNotis config istate) (updateNotisForMe istate)

  atomically $ modifyTVar' istate $
    \istate' -> istate' { stNotiState = notiState }
  createNotiCenter istate config

  unixSignalAdd PRIORITY_HIGH (fromIntegral sigUSR1)
    (showNotiCenter istate notiState)

  runCommand $ configStartupCommand config

  GI.main

main :: IO ()
main = do
--  daemonize main'
  main'


text:
{-# LANGUAGE GADTs #-}

module Benchmark where

import System.Environment

import Ypnos.Examples.Stencils
  (runAvgY, runAvg, runAvgY', raiseToList, runLife, runLife', runId)

import Criterion
import Criterion.Monad
import Criterion.Analysis
import Criterion.Environment
import Criterion.Config
import Statistics.Resampling.Bootstrap
import Data.Vector.Unboxed hiding (map, mapM, foldr, foldr1, (++))
import Prelude hiding (sum, length)
import Control.Monad.Trans

l  = [1,2,3,4,5,6,7,8,9,10]
l2 = [True, False, False, True, False, True]

avgGPU = runAvgY l
avgCPU = runAvgY' l
avgA = raiseToList runAvg l

lifeCPU = runLife' l2
lifeGPU = runLife  l2

idGPU = runId l

stenBench :: ((Int,Int) -> b) -> [Benchmark]
stenBench f = [ bench "10x10" $ whnf f (10,10)
              , bench "100x10" $ whnf f (100,10)
              , bench "100x100" $ whnf f (100,100)
              --, bench "1000x1000" $ whnf f (1000,1000)
              ]

--main = defaultMain [ bgroup "Ypnos" (stenBench runAvgY')
--                   , bgroup "Accelerate" (stenBench runAvgA')
--                   ]
type Result = (Double, Double, Double)

runB :: Fun -> Int -> IO Result
runB (Fun f) x = let v = do env <- measureEnvironment
                            l <- runBenchmark env (whnf f (x,x))
                            cfg <- getConfig
                            let conf = fromLJ cfgConfInterval cfg
                            let res  = fromLJ cfgResamples cfg
                            s <- liftIO $ analyseSample conf l res
                            let m = anMean s
                            return (estPoint m,
                                    estLowerBound m,
                                    estUpperBound m)

         in  withConfig defaultConfig v

makeSet :: (Int -> IO Result) -> [Int] -> IO [[String]]
makeSet f range = let tup x = do (y, l, u)  <- f x
                                 return [show x, show y, show l, show u]
                  in  mapM tup range

insert i x y = x ++ i ++ y

makeLine [] = ""
makeLine xs = foldr1 (insert ", ") xs

printCSV :: [[String]] -> String
printCSV lls = foldr (insert "\n") "" (map makeLine lls)

type FunS = String
type ImplS = String

data Fun where
  Fun :: ((Int,Int) -> b) -> Fun

getFun :: FunS -> ImplS -> Fun
getFun "avg" "gpu" = Fun avgGPU
getFun "avg" "cpu" = Fun avgCPU
getFun "life" "gpu" = Fun lifeGPU
getFun "life" "cpu" = Fun lifeCPU
getFun "id" "gpu" = Fun idGPU

main = do [function, impl, begin, step, end, filename] <- getArgs
          let fun = getFun function impl
          let b = read begin :: Int
          let s = read step :: Int
          let e = read end :: Int

          a <- makeSet (runB fun) ([b, s.. e])
          let ca = printCSV a
          writeFile filename ca

max_stars_repo_name: dorchard/ypnos

max_issues_repo_name: dorchard/ypnos

max_forks_repo_name: dorchard/ypnos


text:
module PlotUtils where

import Data.Complex ( Complex((:+)) )
import Graphics.GD ( newImage, rgb, savePngFile, setPixel, Color, Point, Size )
import Fractals ( iteratingFractal, Formula , flipForJulia)

type ViewPoint = (Float, Float) -- note that while gd's Point will refer to pixel points, ViewPoint will refer to points in the view port, on the complex plane
type ViewPort = (ViewPoint, ViewPoint) -- the exact window frame on the complex plane we want to observe
type JuliaMode = Bool

gridBy :: Size -> [(Int, Int)]
gridBy (width, height) = [(x, y) | x <- [0 .. width - 1], y <- [0 .. height - 1]]

colorBy :: Int -> Color
colorBy x = rgb (x*3+60) (x*2+30) (x^2+ 60) -- change by color ratio wanted

toColor :: ViewPoint -> Formula -> JuliaMode ->  Color
toColor (vx, vy) formula juliaMode = colorBy $ iteratingFractal formula zc 0
  where zc = flipForJulia (0 :+ 0, vx :+ vy) juliaMode

-- pointing from the outside world to a specific point on the complex plane outputing its x,y as a ViewPoint
pixelToViewPoint :: Point -> Size -> ViewPort -> ViewPoint
pixelToViewPoint (pixX, pixY) (width, height) ((minXvp, minYvp), (maxXvp, maxYvp)) =
  (xVP, yVP) where

    xscale = fromIntegral pixX / fromIntegral width
    yscale = fromIntegral pixY / fromIntegral height

    xVP = (maxXvp - minXvp) * xscale + minXvp
    yVP = (maxYvp - minYvp) * yscale + minYvp

plot :: ViewPort -> Size -> Formula -> JuliaMode -> String -> IO ()
plot viewPort size formula juliaMode fileName = do
  img <- newImage size

  mapM_ (\point ->
          let
            viewPoint = pixelToViewPoint point size viewPort
            color = toColor viewPoint formula juliaMode
          in
            setPixel point color img
        ) (gridBy size)

  savePngFile fileName img


text:
{-# LANGUAGE FlexibleInstances #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Data.AEq
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- A type class for approximate and exact equalilty comparisons and instances
-- for common data types.
module Data.AEq (
    AEq(..),
    ) where

import Foreign
import Foreign.C.Types
import Data.Complex
import Numeric.IEEE

infix 4 ===, ~==

-- | Types with approximate and exact equality comparisons.
class Eq a => AEq a where
    -- | An exact equality comparison.
    --
    -- For real 'IEEE' types, two values are equivalent in the
    -- following cases:
    --
    --   * both values are @+0@;
    --
    --   * both values are @-0@;
    --
    --   * both values are nonzero and equal to each other
    --     (according to '==');
    --
    --   * both values are @NaN@ with the same payload and sign.
    --
    -- For complex 'IEEE' types, two values are equivalent if their
    -- real and imaginary parts are equivalent.
    --
    (===) :: a -> a -> Bool
    (===) = (==)
    {-# INLINE (===) #-}

    -- | An approximate equality comparison operator.
    --
    -- For real 'IEEE' types, two values are approximately equal in the
    -- following cases:
    --
    --   * at least half of their significand bits agree;
    --
    --   * both values are less than 'epsilon';
    --
    --   * both values are @NaN@.
    --
    -- For complex 'IEEE' types, two values are approximately equal in the
    -- followiing cases:
    --
    --   * their magnitudes are approximately equal and the angle between
    --     them is less than @32*'epsilon'@;
    --
    --   * both magnitudes are less than 'epsilon';
    --
    --   * both have a @NaN@ real or imaginary part.
    --
    -- Admitedly, the @32@ is a bit of a hack.  Future versions of the
    -- library may switch to a more principled test of the angle.
    --
    (~==) :: a -> a -> Bool
    (~==) = (==)
    {-# INLINE (~==) #-}

approxEqIEEE :: (IEEE a) => a -> a -> Bool
approxEqIEEE x y =
    ( sameSignificandBits x y >= d
    || (abs x < epsilon && abs y < epsilon)
    || (isNaN x && isNaN y)
    )
  where
    d = (floatDigits x + 1) `div` 2
{-# INLINE approxEqIEEE #-}

identicalComplexIEEE :: (IEEE a) => Complex a -> Complex a -> Bool
identicalComplexIEEE (x1 :+ y1) (x2 :+ y2) =
    (identicalIEEE x1 x2) && (identicalIEEE y1 y2)
{-# INLINE identicalComplexIEEE #-}

approxEqComplexIEEE :: (IEEE a) => Complex a -> Complex a -> Bool
approxEqComplexIEEE z1 z2 = let
    (r1,c1) = polar z1
    (r2,c2) = polar z2
    angle = abs (c1 - c2)
    in ( ( approxEqIEEE r1 r2
         && (angle < 32*epsilon || angle > 2*(pi - 16*epsilon) || isNaN angle)
         )
       || (r1 < epsilon && r2 < epsilon)
       )
{-# INLINE approxEqComplexIEEE #-}

instance AEq Float where
    (===) = identicalIEEE
    {-# INLINE (===) #-}
    (~==) = approxEqIEEE
    {-# INLINE (~==) #-}

instance AEq Double where
    (===) = identicalIEEE
    {-# INLINE (===) #-}
    (~==) = approxEqIEEE
    {-# INLINE (~==) #-}

instance AEq (Complex Float) where
    (===) = identicalComplexIEEE
    {-# INLINE (===) #-}
    (~==) = approxEqComplexIEEE
    {-# INLINE (~==) #-}

instance AEq (Complex Double) where
    (===) = identicalComplexIEEE
    {-# INLINE (===) #-}
    (~==) = approxEqComplexIEEE
    {-# INLINE (~==) #-}

instance AEq CFloat where
    (===) = identicalIEEE
    {-# INLINE (===) #-}
    (~==) = approxEqIEEE
    {-# INLINE (~==) #-}

instance AEq CDouble where
    (===) = identicalIEEE
    {-# INLINE (===) #-}
    (~==) = approxEqIEEE
    {-# INLINE (~==) #-}

instance AEq (Complex CFloat) where
    (===) = identicalComplexIEEE
    {-# INLINE (===) #-}
    (~==) = approxEqComplexIEEE
    {-# INLINE (~==) #-}

instance AEq (Complex CDouble) where
    (===) = identicalComplexIEEE
    {-# INLINE (===) #-}
    (~==) = approxEqComplexIEEE
    {-# INLINE (~==) #-}

instance AEq Bool
instance AEq Char
instance AEq Int
instance AEq Int8
instance AEq Int16
instance AEq Int32
instance AEq Int64
instance AEq Integer
instance AEq Ordering
instance AEq Word
instance AEq Word8
instance AEq Word16
instance AEq Word32
instance AEq Word64
instance AEq ()
instance AEq WordPtr
instance AEq IntPtr
instance AEq (StablePtr a)
instance AEq (Ptr a)
instance AEq (FunPtr a)
instance AEq (ForeignPtr a)
instance AEq CChar
instance AEq CSChar
instance AEq CUChar
instance AEq CShort
instance AEq CUShort
instance AEq CInt
instance AEq CUInt
instance AEq CLong
instance AEq CULong
instance AEq CPtrdiff
instance AEq CSize
instance AEq CWchar
instance AEq CSigAtomic
instance AEq CLLong
instance AEq CULLong
instance AEq CIntPtr
instance AEq CUIntPtr
instance AEq CIntMax
instance AEq CUIntMax
instance AEq CClock
instance AEq CTime

eqListsWith :: (a -> a -> Bool) -> [a] -> [a] -> Bool
eqListsWith f (x:xs) (y:ys) = f x y && eqListsWith f xs ys
eqListsWith _ [] [] = True
eqListsWith _ _  _  = False
{-# INLINE eqListsWith #-}

instance (AEq a) => AEq [a] where
    (===) = eqListsWith (===)
    {-# INLINE (===) #-}
    (~==) = eqListsWith (~==)
    {-# INLINE (~==) #-}

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

    (~==) Nothing  Nothing  = True
    (~==) (Just x) (Just y) = (~==) x y
    (~==) _ _ = False
    {-# INLINE (~==) #-}

instance (AEq a, AEq b) => AEq (Either a b) where
    (===) (Left a1)  (Left a2)  = (===) a1 a2
    (===) (Right b1) (Right b2) = (===) b1 b2
    (===) _ _ = False
    {-# INLINE (===) #-}

    (~==) (Left a1)  (Left a2)  = (~==) a1 a2
    (~==) (Right b1) (Right b2) = (~==) b1 b2
    (~==) _ _ = False
    {-# INLINE (~==) #-}

instance (AEq a, AEq b) => AEq (a,b) where
    (===) (a1,b1) (a2,b2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1) (a2,b2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c) => AEq (a,b,c) where
    (===) (a1,b1,c1) (a2,b2,c2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1) (a2,b2,c2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d) => AEq (a,b,c,d) where
    (===) (a1,b1,c1,d1) (a2,b2,c2,d2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1) (a2,b2,c2,d2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d, AEq e) => AEq (a,b,c,d,e) where
    (===) (a1,b1,c1,d1,e1) (a2,b2,c2,d2,e2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        && ((===) e1 e2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1,e1) (a2,b2,c2,d2,e2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        && ((~==) e1 e2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d, AEq e, AEq f) => AEq (a,b,c,d,e,f) where
    (===) (a1,b1,c1,d1,e1,f1) (a2,b2,c2,d2,e2,f2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        && ((===) e1 e2)
        && ((===) f1 f2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1,e1,f1) (a2,b2,c2,d2,e2,f2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        && ((~==) e1 e2)
        && ((~==) f1 f2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d, AEq e, AEq f, AEq g) => AEq (a,b,c,d,e,f,g) where
    (===) (a1,b1,c1,d1,e1,f1,g1) (a2,b2,c2,d2,e2,f2,g2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        && ((===) e1 e2)
        && ((===) f1 f2)
        && ((===) g1 g2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1,e1,f1,g1) (a2,b2,c2,d2,e2,f2,g2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        && ((~==) e1 e2)
        && ((~==) f1 f2)
        && ((~==) g1 g2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d, AEq e, AEq f, AEq g, AEq h) => AEq (a,b,c,d,e,f,g,h) where
    (===) (a1,b1,c1,d1,e1,f1,g1,h1) (a2,b2,c2,d2,e2,f2,g2,h2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        && ((===) e1 e2)
        && ((===) f1 f2)
        && ((===) g1 g2)
        && ((===) h1 h2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1,e1,f1,g1,h1) (a2,b2,c2,d2,e2,f2,g2,h2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        && ((~==) e1 e2)
        && ((~==) f1 f2)
        && ((~==) g1 g2)
        && ((~==) h1 h2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d, AEq e, AEq f, AEq g, AEq h, AEq i) => AEq (a,b,c,d,e,f,g,h,i) where
    (===) (a1,b1,c1,d1,e1,f1,g1,h1,i1) (a2,b2,c2,d2,e2,f2,g2,h2,i2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        && ((===) e1 e2)
        && ((===) f1 f2)
        && ((===) g1 g2)
        && ((===) h1 h2)
        && ((===) i1 i2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1,e1,f1,g1,h1,i1) (a2,b2,c2,d2,e2,f2,g2,h2,i2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        && ((~==) e1 e2)
        && ((~==) f1 f2)
        && ((~==) g1 g2)
        && ((~==) h1 h2)
        && ((~==) i1 i2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d, AEq e, AEq f, AEq g, AEq h, AEq i, AEq j) => AEq (a,b,c,d,e,f,g,h,i,j) where
    (===) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        && ((===) e1 e2)
        && ((===) f1 f2)
        && ((===) g1 g2)
        && ((===) h1 h2)
        && ((===) i1 i2)
        && ((===) j1 j2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        && ((~==) e1 e2)
        && ((~==) f1 f2)
        && ((~==) g1 g2)
        && ((~==) h1 h2)
        && ((~==) i1 i2)
        && ((~==) j1 j2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d, AEq e, AEq f, AEq g, AEq h, AEq i, AEq j, AEq k) => AEq (a,b,c,d,e,f,g,h,i,j,k) where
    (===) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1,k1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2,k2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        && ((===) e1 e2)
        && ((===) f1 f2)
        && ((===) g1 g2)
        && ((===) h1 h2)
        && ((===) i1 i2)
        && ((===) j1 j2)
        && ((===) k1 k2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1,k1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2,k2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        && ((~==) e1 e2)
        && ((~==) f1 f2)
        && ((~==) g1 g2)
        && ((~==) h1 h2)
        && ((~==) i1 i2)
        && ((~==) j1 j2)
        && ((~==) k1 k2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d, AEq e, AEq f, AEq g, AEq h, AEq i, AEq j, AEq k, AEq l) => AEq (a,b,c,d,e,f,g,h,i,j,k,l) where
    (===) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1,k1,l1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2,k2,l2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        && ((===) e1 e2)
        && ((===) f1 f2)
        && ((===) g1 g2)
        && ((===) h1 h2)
        && ((===) i1 i2)
        && ((===) j1 j2)
        && ((===) k1 k2)
        && ((===) l1 l2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1,k1,l1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2,k2,l2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        && ((~==) e1 e2)
        && ((~==) f1 f2)
        && ((~==) g1 g2)
        && ((~==) h1 h2)
        && ((~==) i1 i2)
        && ((~==) j1 j2)
        && ((~==) k1 k2)
        && ((~==) l1 l2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d, AEq e, AEq f, AEq g, AEq h, AEq i, AEq j, AEq k, AEq l, AEq m) => AEq (a,b,c,d,e,f,g,h,i,j,k,l,m) where
    (===) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1,k1,l1,m1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2,k2,l2,m2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        && ((===) e1 e2)
        && ((===) f1 f2)
        && ((===) g1 g2)
        && ((===) h1 h2)
        && ((===) i1 i2)
        && ((===) j1 j2)
        && ((===) k1 k2)
        && ((===) l1 l2)
        && ((===) m1 m2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1,k1,l1,m1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2,k2,l2,m2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        && ((~==) e1 e2)
        && ((~==) f1 f2)
        && ((~==) g1 g2)
        && ((~==) h1 h2)
        && ((~==) i1 i2)
        && ((~==) j1 j2)
        && ((~==) k1 k2)
        && ((~==) l1 l2)
        && ((~==) m1 m2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d, AEq e, AEq f, AEq g, AEq h, AEq i, AEq j, AEq k, AEq l, AEq m, AEq n) => AEq (a,b,c,d,e,f,g,h,i,j,k,l,m,n) where
    (===) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1,k1,l1,m1,n1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2,k2,l2,m2,n2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        && ((===) e1 e2)
        && ((===) f1 f2)
        && ((===) g1 g2)
        && ((===) h1 h2)
        && ((===) i1 i2)
        && ((===) j1 j2)
        && ((===) k1 k2)
        && ((===) l1 l2)
        && ((===) m1 m2)
        && ((===) n1 n2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1,k1,l1,m1,n1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2,k2,l2,m2,n2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        && ((~==) e1 e2)
        && ((~==) f1 f2)
        && ((~==) g1 g2)
        && ((~==) h1 h2)
        && ((~==) i1 i2)
        && ((~==) j1 j2)
        && ((~==) k1 k2)
        && ((~==) l1 l2)
        && ((~==) m1 m2)
        && ((~==) n1 n2)
        )
    {-# INLINE (~==) #-}

instance (AEq a, AEq b, AEq c, AEq d, AEq e, AEq f, AEq g, AEq h, AEq i, AEq j, AEq k, AEq l, AEq m, AEq n, AEq o) => AEq (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o) where
    (===) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1,k1,l1,m1,n1,o1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2,k2,l2,m2,n2,o2) =
        (  ((===) a1 a2)
        && ((===) b1 b2)
        && ((===) c1 c2)
        && ((===) d1 d2)
        && ((===) e1 e2)
        && ((===) f1 f2)
        && ((===) g1 g2)
        && ((===) h1 h2)
        && ((===) i1 i2)
        && ((===) j1 j2)
        && ((===) k1 k2)
        && ((===) l1 l2)
        && ((===) m1 m2)
        && ((===) n1 n2)
        && ((===) o1 o2)
        )
    {-# INLINE (===) #-}

    (~==) (a1,b1,c1,d1,e1,f1,g1,h1,i1,j1,k1,l1,m1,n1,o1) (a2,b2,c2,d2,e2,f2,g2,h2,i2,j2,k2,l2,m2,n2,o2) =
        (  ((~==) a1 a2)
        && ((~==) b1 b2)
        && ((~==) c1 c2)
        && ((~==) d1 d2)
        && ((~==) e1 e2)
        && ((~==) f1 f2)
        && ((~==) g1 g2)
        && ((~==) h1 h2)
        && ((~==) i1 i2)
        && ((~==) j1 j2)
        && ((~==) k1 k2)
        && ((~==) l1 l2)
        && ((~==) m1 m2)
        && ((~==) n1 n2)
        && ((~==) o1 o2)
        )
    {-# INLINE (~==) #-}


text:
{-# LANGUAGE FlexibleContexts #-}
-- |
-- Module    : Statistics.Test.KolmogorovSmirnov
-- Copyright : (c) 2011 Aleksey Khudyakov
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Kolmogov-Smirnov tests are non-parametric tests for assessing
-- whether given sample could be described by distribution or whether
-- two samples have the same distribution. It's only applicable to
-- continuous distributions.
module Statistics.Test.KolmogorovSmirnov (
    -- * Kolmogorov-Smirnov test
    kolmogorovSmirnovTest
  , kolmogorovSmirnovTestCdf
  , kolmogorovSmirnovTest2
    -- * Evaluate statistics
  , kolmogorovSmirnovCdfD
  , kolmogorovSmirnovD
  , kolmogorovSmirnov2D
    -- * Probablities
  , kolmogorovSmirnovProbability
    -- * References
    -- $references
  , module Statistics.Test.Types
  ) where

import Control.Monad (when)
import Prelude hiding (exponent, sum)
import Statistics.Distribution (Distribution(..))
import Statistics.Function (gsort, unsafeModify)
import Statistics.Matrix (center, for, fromVector)
import qualified Statistics.Matrix as Mat
import Statistics.Test.Types
import Statistics.Types (mkPValue)
import qualified Data.Vector          as V
import qualified Data.Vector.Storable as S
import qualified Data.Vector.Unboxed  as U
import qualified Data.Vector.Generic  as G
import           Data.Vector.Generic    ((!))
import qualified Data.Vector.Unboxed.Mutable as M

----------------------------------------------------------------
-- Test
----------------------------------------------------------------

-- | Check that sample could be described by distribution. Returns
--   @Nothing@ is sample is empty
--
--   This test uses Marsaglia-Tsang-Wang exact algorithm for
--   calculation of p-value.
kolmogorovSmirnovTest :: (Distribution d, G.Vector v Double)
                      => d        -- ^ Distribution
                      -> v Double -- ^ Data sample
                      -> Maybe (Test ())
{-# INLINE kolmogorovSmirnovTest #-}
kolmogorovSmirnovTest d
  = kolmogorovSmirnovTestCdf (cumulative d)

-- | Variant of 'kolmogorovSmirnovTest' which uses CDF in form of
--   function.
kolmogorovSmirnovTestCdf :: (G.Vector v Double)
                         => (Double -> Double) -- ^ CDF of distribution
                         -> v Double           -- ^ Data sample
                         -> Maybe (Test ())
{-# INLINE kolmogorovSmirnovTestCdf #-}
kolmogorovSmirnovTestCdf cdf sample
  | G.null sample = Nothing
  | otherwise     = Just Test
      { testSignificance = mkPValue $ 1 - prob
      , testStatistics   = d
      , testDistribution = ()
      }
  where
    d    = kolmogorovSmirnovCdfD cdf sample
    prob = kolmogorovSmirnovProbability (G.length sample) d

-- | Two sample Kolmogorov-Smirnov test. It tests whether two data
--   samples could be described by the same distribution without
--   making any assumptions about it. If either of samples is empty
--   returns Nothing.
--
--   This test uses approximate formula for computing p-value.
kolmogorovSmirnovTest2 :: (G.Vector v Double)
                       => v Double -- ^ Sample 1
                       -> v Double -- ^ Sample 2
                       -> Maybe (Test ())
kolmogorovSmirnovTest2 xs1 xs2
  | G.null xs1 || G.null xs2 = Nothing
  | otherwise                = Just Test
      { testSignificance = mkPValue $ 1 - prob d
      , testStatistics   = d
      , testDistribution = ()
      }
  where
    d    = kolmogorovSmirnov2D xs1 xs2
         * (en + 0.12 + 0.11/en)
    -- Effective number of data points
    n1   = fromIntegral (G.length xs1)
    n2   = fromIntegral (G.length xs2)
    en   = sqrt $ n1 * n2 / (n1 + n2)
    --
    prob z
      | z <  0    = error "kolmogorovSmirnov2D: internal error"
      | z == 0    = 0
      | z <  1.18 = let y = exp( -1.23370055013616983 / (z*z) )
                    in  2.25675833419102515 * sqrt( -log y ) * (y + y**9 + y**25 + y**49)
      | otherwise = let x = exp(-2 * z * z)
                    in  1 - 2*(x - x**4 + x**9)
{-# INLINABLE  kolmogorovSmirnovTest2 #-}
{-# SPECIALIZE kolmogorovSmirnovTest2 :: U.Vector Double -> U.Vector Double -> Maybe (Test ()) #-}
{-# SPECIALIZE kolmogorovSmirnovTest2 :: V.Vector Double -> V.Vector Double -> Maybe (Test ()) #-}
{-# SPECIALIZE kolmogorovSmirnovTest2 :: S.Vector Double -> S.Vector Double -> Maybe (Test ()) #-}
-- FIXME: Find source for approximation for D

----------------------------------------------------------------
-- Kolmogorov's statistic
----------------------------------------------------------------

-- | Calculate Kolmogorov's statistic /D/ for given cumulative
--   distribution function (CDF) and data sample. If sample is empty
--   returns 0.
kolmogorovSmirnovCdfD :: G.Vector v Double
                      => (Double -> Double) -- ^ CDF function
                      -> v Double           -- ^ Sample
                      -> Double
kolmogorovSmirnovCdfD cdf sample
  | G.null sample = 0
  | otherwise     = G.maximum
                  $ G.zipWith3 (\p a b -> abs (p-a) `max` abs (p-b))
                    ps steps (G.tail steps)
  where
    xs = gsort sample
    n  = G.length xs
    --
    ps    = G.map cdf xs
    steps = G.map (/ fromIntegral n)
          $ G.generate (n+1) fromIntegral
{-# INLINABLE  kolmogorovSmirnovCdfD #-}
{-# SPECIALIZE kolmogorovSmirnovCdfD :: (Double -> Double) -> U.Vector Double -> Double #-}
{-# SPECIALIZE kolmogorovSmirnovCdfD :: (Double -> Double) -> V.Vector Double -> Double #-}
{-# SPECIALIZE kolmogorovSmirnovCdfD :: (Double -> Double) -> S.Vector Double -> Double #-}

-- | Calculate Kolmogorov's statistic /D/ for given cumulative
--   distribution function (CDF) and data sample. If sample is empty
--   returns 0.
kolmogorovSmirnovD :: (Distribution d, G.Vector v Double)
                   => d         -- ^ Distribution
                   -> v Double  -- ^ Sample
                   -> Double
kolmogorovSmirnovD d = kolmogorovSmirnovCdfD (cumulative d)
{-# INLINE kolmogorovSmirnovD #-}

-- | Calculate Kolmogorov's statistic /D/ for two data samples. If
--   either of samples is empty returns 0.
kolmogorovSmirnov2D :: (G.Vector v Double)
                    => v Double   -- ^ First sample
                    -> v Double   -- ^ Second sample
                    -> Double
kolmogorovSmirnov2D sample1 sample2
  | G.null sample1 || G.null sample2 = 0
  | otherwise                        = worker 0 0 0
  where
    xs1 = gsort sample1
    xs2 = gsort sample2
    n1  = G.length xs1
    n2  = G.length xs2
    en1 = fromIntegral n1
    en2 = fromIntegral n2
    -- Find new index
    skip x i xs = go (i+1)
      where go n | n >= G.length xs = n
                 | xs ! n == x      = go (n+1)
                 | otherwise        = n
    -- Main loop
    worker d i1 i2
      | i1 >= n1 || i2 >= n2 = d
      | otherwise            = worker d' i1' i2'
      where
        d1  = xs1 ! i1
        d2  = xs2 ! i2
        i1' | d1 <= d2  = skip d1 i1 xs1
            | otherwise = i1
        i2' | d2 <= d1  = skip d2 i2 xs2
            | otherwise = i2
        d'  = max d (abs $ fromIntegral i1' / en1 - fromIntegral i2' / en2)
{-# INLINABLE  kolmogorovSmirnov2D #-}
{-# SPECIALIZE kolmogorovSmirnov2D :: U.Vector Double -> U.Vector Double -> Double #-}
{-# SPECIALIZE kolmogorovSmirnov2D :: V.Vector Double -> V.Vector Double -> Double #-}
{-# SPECIALIZE kolmogorovSmirnov2D :: S.Vector Double -> S.Vector Double -> Double #-}

-- | Calculate cumulative probability function for Kolmogorov's
--   distribution with /n/ parameters or probability of getting value
--   smaller than /d/ with n-elements sample.
--
--   It uses algorithm by Marsgalia et. al. and provide at least
--   7-digit accuracy.
kolmogorovSmirnovProbability :: Int    -- ^ Size of the sample
                             -> Double -- ^ D value
                             -> Double
kolmogorovSmirnovProbability n d
  -- Avoid potentially lengthy calculations for large N and D > 0.999
  | s > 7.24 || (s > 3.76 && n > 99) = 1 - 2 * exp( -(2.000071 + 0.331 / sqrt n' + 1.409 / n') * s)
  -- Exact computation
  | otherwise = fini $ KSMatrix 0 matrix `power` n
  where
    s  = n' * d * d
    n' = fromIntegral n

    size = 2*k - 1
    k    = floor (n' * d) + 1
    h    = fromIntegral k - n' * d
    -- Calculate initial matrix
    matrix =
      let m = U.create $ do
            mat <- M.new (size*size)
            -- Fill matrix with 0 and 1s
            for 0 size $ \row ->
              for 0 size $ \col -> do
                let val | row + 1 >= col = 1
                        | otherwise      = 0 :: Double
                M.write mat (row * size + col) val
            -- Correct left column/bottom row
            for 0 size $ \i -> do
              let delta = h ^^ (i + 1)
              unsafeModify mat (i    * size)         (subtract delta)
              unsafeModify mat (size * size - 1 - i) (subtract delta)
            -- Correct corner element if needed
            when (2*h > 1) $ do
              unsafeModify mat ((size - 1) * size) (+ ((2*h - 1) ^ size))
            -- Divide diagonals by factorial
            let divide g num
                  | num == size = return ()
                  | otherwise   = do for num size $ \i ->
                                       unsafeModify mat (i * (size + 1) - num) (/ g)
                                     divide (g * fromIntegral (num+2)) (num+1)
            divide 2 1
            return mat
      in fromVector size size m
    -- Last calculation
    fini (KSMatrix e m) = loop 1 (center m) e
      where
        loop i ss eQ
          | i  > n       = ss * 10 ^^ eQ
          | ss' < 1e-140 = loop (i+1) (ss' * 1e140) (eQ - 140)
          | otherwise    = loop (i+1)  ss'           eQ
          where ss' = ss * fromIntegral i / fromIntegral n

data KSMatrix = KSMatrix Int Mat.Matrix

multiply :: KSMatrix -> KSMatrix -> KSMatrix
multiply (KSMatrix e1 m1) (KSMatrix e2 m2) = KSMatrix (e1+e2) (Mat.multiply m1 m2)

power :: KSMatrix -> Int -> KSMatrix
power mat 1 = mat
power mat n = avoidOverflow res
  where
    mat2 = power mat (n `quot` 2)
    pow  = multiply mat2 mat2
    res | odd n     = multiply pow mat
        | otherwise = pow

avoidOverflow :: KSMatrix -> KSMatrix
avoidOverflow ksm@(KSMatrix e m)
  | center m > 1e140 = KSMatrix (e + 140) (Mat.map (* 1e-140) m)
  | otherwise        = ksm

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

-- $references
--
-- * G. Marsaglia, W. W. Tsang, J. Wang (2003) Evaluating Kolmogorov's
--   distribution, Journal of Statistical Software, American
--   Statistical Association, vol. 8(i18).


text:
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Signal
-- Copyright   :  (c) Alexander Vivian Hugh McPhail 2010, 2014, 2015, 2016
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  uses FFI
--
-- Signal processing functions
--
-----------------------------------------------------------------------------

module Numeric.Signal (
                       S.Convolvable(..),
                       S.Filterable(),
                       -- * Filtering
                       hamming,
                       pwelch,
                       fir,standard_fir,broadband_fir,
                       freqzF,freqzN,
                       filter,broadband_filter,
                       -- * Analytic Signal
                       analytic_signal,analytic_power,analytic_phase,
                       unwrap,
                       -- * Statistics
                       cross_covariance,cross_correlation,cross_spectrum,
                       auto_covariance,auto_correlation,
                       -- * Preprocessing
                       detrend,
                       resize,
                       downsample,
                       deriv,
                       cumulative_sum
                ) where

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

import qualified Numeric.Signal.Internal as S

import Numeric.GSL.Fitting.Linear

import Data.Complex
import Foreign.Storable()

--import Data.Function.Unicode

import qualified Data.List as L

--import Data.Packed.Vector
--import Data.Packed(Container(..))
import Numeric.LinearAlgebra

import qualified Data.Vector.Generic as GV

import qualified Numeric.GSL.Fourier as F

import Prelude hiding(filter)

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

-- | filters the signal
filter :: (S.Filterable a)
         => Vector a   -- ^ zero coefficients
       -> Vector a     -- ^ pole coefficients
       -> Int   -- ^ sampling rate
       -> Vector a     -- ^ input signal
       -> Vector a     -- ^ output signal
filter b a s v = let len = size v
                     w = min s len
                     start = (negate . fromList . reverse . toList . subVector 0 w) v
                     finish = (negate . fromList . reverse . toList . subVector (len-w) w) v
                     v' = vjoin [start,v,finish]
                 in subVector s len $ S.filter_ b a v'

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

-- | Welch (1967) power spectrum density using periodogram/FFT method
pwelch :: Int            -- ^ sampling rate
       -> Int            -- ^ window size
       -> Vector Double  -- ^ input signal
       -> (Vector Double,Vector Double)  -- ^ (frequency index,power density)
pwelch s w v = let w' = max s w -- make window at least sampling rate
                   r  = S.pwelch w' v
                   sd = (fromIntegral s)/2
                   -- scale for sampling rate
                   r' = scale (recip sd) r
                   f  = linspace ((w `div` 2) + 1) (0,sd)
               in (f,r')

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

-- | a broadband FIR
broadband_fir :: (S.Filterable a, Double ~ DoubleOf a, Convert (Complex a)) =>
                Int           -- ^ sampling rate
              -> (Int,Int)     -- ^ (lower,upper) frequency cutoff
              -> Vector a -- ^ filter coefficients
broadband_fir s (l,h) = let o = 501
                            ny = (fromIntegral s) / 2.0
                            fl = (fromIntegral l) / ny
                            fh = (fromIntegral h) / ny
                            f = [0, fl*0.95, fl, fh, fh*1.05, 1]
                            m = [0,0,1,1,0,0]
                            be = zip f m
                        in standard_fir o be

-- | a broadband filter
broadband_filter :: (S.Filterable a, Double ~ DoubleOf a)
                   => Int        -- ^ sampling rate
                 -> (Int,Int)    -- ^ (lower,upper) frequency cutoff
                 -> Vector a            -- ^ input signal
                 -> Vector a            -- ^ output signal
broadband_filter s f v = let b = S.fromDouble $ broadband_fir s f
                         in filter b (scalar 1.0) s v

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

-- | standard FIR filter
-- |   FIR filter with grid a power of 2 greater than the order, ramp = grid/16, hamming window
standard_fir :: (S.Filterable a, Double ~ DoubleOf a, Convert (Complex a)) =>
               Int -> [(a,a)] -> Vector a
standard_fir o be = let grid  = calc_grid o
                        trans_ = grid `div` 16
                    in fir o be grid trans_ $ S.hamming_ (o+1)

calc_grid :: Int -> Int
calc_grid o = let next_power = ceiling (((log $ fromIntegral o) :: Double) / (log 2.0)) :: Int
              in floor $ 2.0 ** ((fromIntegral next_power) :: Double)

-- | produce an FIR filter
fir :: (S.Filterable a
      , Convert (Complex a), Double ~ DoubleOf a) =>
      Int               -- ^ order (one less than the length of the filter)
    -> [(a,a)] -- ^ band edge frequency, nondecreasing, [0, f1, ..., f(n-1), 1]
                        -- ^ band edge magnitude
    -> Int               -- ^ grid spacing
    -> Int               -- ^ transition width
    -> Vector a     -- ^ smoothing window (size is order + 1)
    -> Vector a     -- ^ the filter coefficients
fir o be gn tn w = let mid = o `div` 2
                       (f,m) = unzip be
                       f' = diff (((fromIntegral gn))/((fromIntegral tn))/2.0) f
                       m' = interpolate f m f'
                       grid = interpolate f' m' $ map (\x -> (fromIntegral x)/(fromIntegral gn)) [0..(gn-1)]
                       grid' = map (\x -> x :+ 0) grid
                       b = S.fromDouble $ fst $ fromComplex $ F.ifft $ double $ fromList $ grid' ++ (reverse (drop 1 grid'))
                       b' = vjoin [subVector ((size b)-mid-1) (mid+1) b, subVector 1 (mid+1) b]
                   in b' * w

floor_zero x
    | x < 0.0   = 0.0
    | otherwise = x

ceil_one x
    | x > 1.0   = 1.0
    | otherwise = x

diff :: S.Filterable a => a -> [a] -> [a]
diff _ []  = []
diff _ [x] = [x]
diff inc (x1:x2:xs)
     | x1 == x2     = (floor_zero $ x1-inc):x1:(ceil_one $ x1+inc):(diff inc (L.filter (/= x2) xs))
     | otherwise    = x1:(diff inc (x2:xs))

interpolate :: S.Filterable a => [a] -> [a] -> [a] -> [a]
interpolate _ _ []      = []
interpolate x y (xp:xs) = if xp == 1.0
                             then ((interpolate'' ((length x)-1) x y xp):(interpolate x y xs))
                             else ((interpolate' x y xp):(interpolate x y xs))

interpolate' :: S.Filterable a => [a] -> [a] -> a -> a
interpolate' x y xp = let Just j = L.findIndex (> xp) x
                      in (interpolate'' j x y xp)

interpolate'' :: S.Filterable a => Int -> [a] -> [a] -> a -> a
interpolate'' j x y xp = let x0 = x !! (j-1)
                             y0 = y !! (j-1)
                             x1 = x !! j
                             y1 = y !! j
                         in y0 + (xp - x0) * ((y1 - y0)/(x1-x0))

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

-- | determine the frequency response of a filter, given a vector of frequencies
freqzF :: (S.Filterable a, Double ~ DoubleOf a, S.Filterable (DoubleOf a)) =>
         Vector a     -- ^ zero coefficients
       -> Vector a       -- ^ pole coefficients
       -> Int     -- ^ sampling rate
       -> Vector a       -- ^ frequencies
       -> Vector a       -- ^ frequency response
freqzF b a s f = S.freqz b a ((2*pi/(fromIntegral s)) * f)

-- | determine the frequency response of a filter, given a number of points and sampling rate
freqzN :: (S.Filterable a, Double ~ DoubleOf a) =>
         Vector a     -- ^ zero coefficients
       -> Vector a       -- ^ pole coefficients
       -> Int     -- ^ sampling rate
       -> Int     -- ^ number of points
       -> (Vector a,Vector a)   -- ^ (frequencies,response)
freqzN b a s n = let w' = linspace n (0,((fromIntegral n)-1)/(fromIntegral (2*n)))
                     r = S.freqz b a ((2*pi)*w')
                     in ((fromIntegral s)*w',r)

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

-- | an analytic signal is the original signal with Hilbert-transformed signal as imaginary component
analytic_signal :: Vector Double -> Vector (Complex Double)
analytic_signal = S.hilbert

-- | the power (amplitude^2 = v * (conj c)) of an analytic signal
analytic_power :: S.Filterable a => Vector (Complex Double) -> Vector a
analytic_power = S.complex_power_

-- | the phase of an analytic signal
analytic_phase :: (S.Filterable a) =>
                 Vector (Complex a) -> Vector a
analytic_phase = (uncurry arctan2) . fromComplex

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

-- | remove a linear trend from data
detrend :: Int             -- ^ window size
        -> Vector Double   -- ^ data to be detrended
        -> Vector Double   -- ^ detrended data
detrend w v = let windows = size v `div` w
                  re = size v - (windows * w)
                  re' = if re == 0 then [] else [re]
                  ws = takesV ((replicate windows w) ++ re') v
                  ds = map detrend' ws
                  windows' = (size v - (w `div` 2)) `div` w
                  ws' = takesV (((w `div` 2):(replicate windows' w)) ++ [size v - (w `div` 2) - (windows' * w)]) v
                  ds' = map detrend' ws'
              in (vjoin ds + vjoin ds') / 2
    where detrend' x = let ln = size x
                           t = linspace ln (1.0,fromIntegral ln)
                           (c0,c1,_,_,_,_) = linear t x
                       in x - (scale c1 t + scalar c0)

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

-- | resize the vector to length n by resampling
resize :: S.Filterable a => Int -> Vector a -> Vector a
resize n v = S.downsample_ (size v `div` n) v

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

-- | cross covariance of two signals
--     the cross correlation is computed by dividing the result
--     by the product of the two standard deviations
cross_covariance :: S.Filterable a =>
                   Int           -- ^ maximum delay
                 -> Vector a -- ^ time series
                 -> Vector a -- ^ time series
                 -> (a,a,Vector a) -- ^ (sd_x,sd_y,cov_xy)
cross_covariance = S.cross_covariance_

-- | cross correlation of two signals
cross_correlation :: S.Filterable a =>
                   Int           -- ^ maximum delay
                 -> Vector a -- ^ time series
                 -> Vector a -- ^ time series
                 -> Vector a -- ^ result
cross_correlation l x y = let (sx,sy,r) = S.cross_covariance_ l x y
                          in GV.map(/ (sx*sy)) r

-- | compute the cross spectrum
cross_spectrum :: (S.Filterable a, Double ~ DoubleOf a) =>
                 Int                      -- ^ maximum delay
               -> Vector a                 -- ^ time series
               -> Vector a                 -- ^ time series
               -> Vector (Complex Double)  -- ^ result
cross_spectrum l x y = (\(_,_,c) -> F.fft (complex $ double c)) (cross_covariance l x y)

-- | auto covariance of two signals
--     the auto correlation is computed by dividing the result
--     by the variance
auto_covariance :: S.Filterable a =>
                   Int           -- ^ maximum delay
                 -> Vector a -- ^ time series
                 -> (a,Vector a) -- ^ (var,cov_xx)
auto_covariance l v = let (sd,_,r) = cross_covariance l v v
                      in (sd*sd,r)

-- | auto correlation of two signals
auto_correlation :: S.Filterable a =>
                   Int           -- ^ maximum delay
                 -> Vector a -- ^ time series
                 -> Vector a -- ^ result
auto_correlation l v = let (var,r) = auto_covariance l v
                          in GV.map(/ var) r

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

-- | coefficients of a Hamming window
hamming :: S.Filterable a =>
          Int           -- ^ length
        -> Vector a -- ^ the Hamming coeffficents
hamming = S.hamming_

-- | resample, take one sample every n samples in the original
downsample :: S.Filterable a => Int -> Vector a -> Vector a
downsample = S.downsample_

-- | the difference between consecutive elements of a vector
deriv :: S.Filterable a => Vector a -> Vector a
deriv = S.deriv_

-- | cumulative sum of a series
cumulative_sum :: S.Filterable a =>
                 Vector a
               -> Vector a
cumulative_sum = S.cumulative_sum_

-- | unwrap the phase of signal (input expected to be within (-pi,pi))
unwrap :: S.Filterable a => Vector a -> Vector a
unwrap = S.unwrap_

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

max_stars_repo_name: amcphail/hsignal

max_issues_repo_name: amcphail/hsignal

max_forks_repo_name: amcphail/hsignal


text:
module Utils where

import Data.Random.Normal
import Numeric.LinearAlgebra
import System.Random
--import Test.QuickCheck

genNormal :: (Random a, Floating a) => [a]
genNormal = let g = mkStdGen 1001 in
  normals g

testGenNormal :: (Random a, Floating a) => [a]
testGenNormal = let g = mkStdGen 1001 in
  take 10 (normals g)

getAverageMat :: Matrix R -> Double
getAverageMat mat =
  let r_f = fromIntegral $ rows mat
      c_f = fromIntegral $ cols mat
      sumMat = sumElements mat in
  sumMat / (r_f * c_f)

max_stars_repo_name: neutronest/vortex

max_issues_repo_name: neutronest/vortex

max_forks_repo_name: neutronest/vortex


text:
{-# LANGUAGE FlexibleInstances #-}

import Debug.Trace
-- import Control.Monad.State.Lazy

import System.Environment (getArgs)
import System.IO (readFile)
import Data.Map.Strict (Map, (!), insert, elems, fromList, toList, findWithDefault, size, empty, member, findMin, findMax, singleton, filter)
import qualified Data.Map.Strict as M
import Data.Set (Set)
import qualified Data.Set as S
--import qualified Data.Array as A
import Data.List (find, intercalate, intersperse, permutations, inits, tails, isPrefixOf)
import Data.List.Split (splitOn)
-- import Data.Complex (Complex((:+)), realPart, imagPart) -- define my own complex

import UI.NCurses
import Data.Char (chr, ord)
-- import Data.Complex (Complex((:+)))
import Data.Maybe (fromJust)
type Instructions = Map Integer Integer

data ComputerState = Ready | Running | Blocked | Done
  deriving (Show, Eq)

data Computer = Computer {
  state :: ComputerState,
  memory :: Map Integer Integer,
  iptr :: Integer, -- instruction pointer
  base :: Integer, -- base offset
  input :: [Integer],
  output :: [Integer]}

computer0 = Computer { state = Ready, memory = empty, iptr = 0, base = 0, input = [], output = [] }

run :: Computer -> Computer -- (Integer, Integer) -> Instructions -> [Integer] -> [Integer]
run c
  | state c == Ready = run $ c {state = Running}
  | state c == Done = c
  | state c == Blocked = if input c == [] then c else run $ c {state = Running}
  | otherwise = -- running
    --traceShow (iptr c, base c, memory c) $
    case instr `mod` 100 of
      1 -> -- add
        run $ c {iptr = i+4, memory = insert (addr 3) (arg 1 + arg 2) instructions}
      2 -> -- multiply
        run $ c {iptr = i+4, memory = insert (addr 3) (arg 1 * arg 2) instructions}
      3 -> -- read input
        if (null $ input c)
        then c {state = Blocked}
        else run $ c {iptr = i+2, memory = insert (addr 1) (head $ input c) instructions, input = tail $ input c}
      4 -> -- output
        --traceShow (arg 1) $
        run $ c {iptr = i+2, output = output c ++ [arg 1]}
      5 -> -- jump-if-true
        run $ c {iptr = if arg 1 == 0 then i+3 else arg 2}
      6 -> -- jump-if-false
        run $ c {iptr = if arg 1 == 0 then arg 2 else i+3}
      7 -> -- less than
        run $ c {iptr = i+4, memory = insert (addr 3) (if arg 1 < arg 2 then 1 else 0) instructions}
      8 -> -- equals
        run $ c {iptr = i+4, memory = insert (addr 3) (if arg 1 == arg 2 then 1 else 0) instructions}
      9 -> -- set relative base
        run $ c {iptr = i+2, base = base c + arg 1}
      99 -> -- halt
        run $ c {state = Done}
      _ -> error "unknown opcode"
  where instructions = memory c
        i = iptr c
        instr = instructions!i
        ii x = findWithDefault 0 x instructions
        arg :: Integer -> Integer
        arg n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ ii $ i+n
          1 -> ii $ i+n
          2 -> ii ((ii $ i+n) + base c)
          _ -> error "bad argument mode"
        addr n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ i+n
          1 -> error "address in mode 1"
          2 -> --trace "address in mode 2" $
            ii (i+n) + base c
          _ -> error $ "address in unknown mode"

-- instance (Show a) => Show (Complex a) where
--   show (C x y) = "C " ++ show x ++ " " ++ show y

-- instance (Eq a) => Eq (Complex a) where
--   (C a ai) == (C b bi) = a == b && ai == bi

-- instance (Eq a, Ord a) => Ord (Complex a) where
--   compare (C a ai) (C b bi) | a < b || a == b && ai < bi = LT
--                             | a > b || a == b && ai > bi = GT
--                             | otherwise = EQ

-- -- Note that we are using a plane where imaginary (y-axis) numbers
-- -- increase downward. I think this means that right and left are
-- -- reversed.

-- data Complex a = C a a

-- -- instance Num (Complex Int) where
-- --   abs c = error "Can't take abs of Complex Int"

-- instance Num a => Num (Complex a) where
--   (C x y) + (C u v) = C (x+u) (y+v)
--   (C x y) * (C u v) = C (x*u-y*v) (x*v+y*u)
--   fromInteger n = C (fromInteger n) 0
--   abs (C x y) = C (abs x) y -- sqrt $ x*x + y*y
--   signum (C x y) = C (signum x) 0
--   negate (C x y) = C (negate x) (negate y)

main = do
  -- [instructionFile] <- getArgs
  instructionStrings <- readFile "21.input.txt" -- instructionFile
  let instructions = fromList . zip [0 ..] $ map read $ splitOn "," instructionStrings

  putStrLn "Part 1"
  let linesToAscii = map (fromIntegral . ord) . unlines
  -- Jumps are 4 spaces. Jump on an island if you can. # = ground, . = hole, _ = don't care
  -- Jump if any of the three spaces in front of you have holes and the fourth space is solid.
  -- @ A B C D
  -- _ . _ _ #
  -- _ _ . _ #
  -- _ _ _ . #
  -- ~(ABC)D -> ~((ABC)+D)
  let ls = [ "OR A J"
           , "AND B J"
           , "AND C J"
           , "NOT J J"
           , "AND D J"
           , "WALK"]

  -- let c = run $ computer0{memory = instructions, input=ls}
  let c = run $ computer0{memory = instructions}
  putStrLn $ map (chr . fromIntegral) $ output c
  putStrLn $ "Running with input: \n" ++ unlines ls
  let c2 = run $ c {input=linesToAscii ls, output=[]}
  let (msg, result) = span (< 128) $ output c2
  putStrLn $ map (chr . fromIntegral) msg ++ concatMap show result

  putStrLn "\n\nPart 2"
  -- Jump as before . . .
  -- @ABCD  EFGHI
  -- _.__#
  -- __._#
  -- ___.#
  -- But check that you are not blocked after the jump.
  --        EFGHI
  --        .__#_
  --        #___#
  --        ##___
  -- ~(ABC)D -> ~((ABC)+D)
  -- (~E)H + EI + EF = (~E)H + E(I+F) = E(F+I) + H

  let ls2 = [
              -- jump if a hole is ahead and the landing is open
              "OR A J"
            , "AND B J"
            , "AND C J"
            , "NOT J J"
            , "AND D J"
              -- and if the future looks ok
            , "OR F T"
            , "OR I T"
            , "AND E T"
            , "OR H T"
            , "AND T J"
            , "RUN"
            ]
  let c3 = run $ computer0{memory = instructions}
  putStrLn $ map (chr . fromIntegral) $ output c3
  putStrLn $ "Running with input: \n" ++ unlines ls2
  let c4 = run $ c3 {input=linesToAscii ls2, output=[]}
  let (msg, result) = span (< 128) $ output c4
  putStrLn $ map (chr . fromIntegral) msg ++ concatMap show result


text:
module Quant.Models (
    CharFunc(..)
) where

import Data.Complex
import Quant.Math.Integration
import Quant.ContingentClaim(OptionType(..))

{- | The 'CharFunc' class defines those
models which have closed-form characteristic
functions.

Minimal complete definition: 'charFunc'.

Still under construction.
-}
class CharFunc a where
    -- | Creates a characteristic function for a model.
    charFunc :: CharFunc a => a -> Double -> Complex Double -> Complex Double

    charFuncOption :: (CharFunc a) => a -> OptionType -> Integrator -> Double
        -> Double -> Double -> Double
    charFuncOption model opt intF strike tmat damp = case opt of
        Put  -> s * q * p1 - k * disc * p2
        Call -> k * disc * (1-p2) - s * q * (1-p1)
            where
                s =
                p1 = 0.5 + intF func1 1E-8 100 / pi
                p2 = 0.5 + intF func2 1E-8 100 / pi

                f v' = realPart $ exp (i*v*k) * leftTerm * rightTerm
                    where
                        v = v' :+ 0
                        damp' = damp :+ 0
                        k = log strike :+ 0
                        i = 0 :+ 1
                        leftTerm = d / (damp' + i * v) / (damp'+i*v+(1:+0))
                        rightTerm = cf $ v - i * (damp' + 1)
                        d = disc yc tmat :+ 0
                        cf x = charFuncMart model fg tmat x

max_stars_repo_name: cdituri/quantfin

max_issues_repo_name: cdituri/quantfin

max_forks_repo_name: cdituri/quantfin


text:
module Augmentation.DGA
    (DGA_Map (DGA_Map)
    ,Augmentation (Aug)
    ,applyDGAMap
    ,compose_maps
    ,Algebra
    ,fromDGAMap
    ) where

import Algebra
import Braid
import Data.Maybe
import Numeric.LinearAlgebra as N
import Data.List
import Debug.Trace

default (Int,Double)

data DGA_Map = DGA_Map [(Char,Algebra)] deriving Eq
instance Show DGA_Map where
    show (DGA_Map l) = "[" ++ (foldr (\(c,e) xs -> [c] ++ "→" ++ show e ++ (if xs /= "]" then "," else "") ++ xs) "]" l)

data Augmentation = Aug StdBraid [(Char,[Vector Z])]
instance Show Augmentation where
    show (Aug b m) = (show b) ++ "\n" ++ (concat $ map (\(c,vs) -> [c] ++ "→" ++ (foldr (\x xs -> show x ++ if xs == "" then "" else "+"++xs) "" vs) ++ " ") m)
instance Eq Augmentation where
    (Aug b1 m1) == (Aug b2 m2) = b1 == b2 && inZ2 m1 == inZ2 m2 && (m1 == m2 || eqh m1 m2)

inZ2 :: [(Char,[Vector Z])] -> [(Char,Z2)]
inZ2 = map (\(c,vs) -> (c,sum $ map (\_ -> fromInteger 1) vs))

isUpperTri :: Matrix Z -> Bool
isUpperTri mat = uptrih mat 0

uptrih :: Matrix Z -> Int -> Bool
uptrih mat i = let (l,w) = size mat
                   cond = and $ map (\j -> (==0) $ mat `atIndex` (j,i)) [i+1..l-1]
                in if i == w then True else cond && (uptrih mat $ i+1)

eqh :: [(Char,[Vector Z])] -> [(Char,[Vector Z])] -> Bool
eqh l1 l2 = maybe False id $ do
            { l1' <- mapM (\c -> lookup c l1) $ map fst l2
            ; l2' <- mapM (\(l,l') -> if length l == length l' then Just l' else Nothing) $ zip l1' $ map snd l2
            ; return True
            ; let dims1 = nub $ map size $ concat l1'
            ; let dims2 = nub $ map size $ concat l2'
            ; dim <- if length dims1 == 1 && length dims2 == 1 && dims1 == dims2 then Just $ head dims1 else Nothing
            ; let m2 = fromRows $ map (fromZ :: Vector Z -> Vector R) $ concat l2'
            ; let n0s = map (\l -> (0,(length $ permutations l)-1)) l1'
            ; let suc mns = if mns == [] then [] else if (fst $ head mns) == (snd $ head mns) then (0,snd $ head mns):(suc $ tail mns) else (1 + (fst $ head mns), snd $ head mns):(tail mns)
            ; let bound (acc,mns) = if mns == [] then acc else bound (((snd $ head mns)+1)*acc,tail mns)
            ; let ubound = bound (1,n0s)
            ; let mat mns = fromRows $ map (fromZ :: Vector Z -> Vector R) $ concat $ zipWith (\l (m,_) -> (permutations l) !! m) l1' mns
            ; let check m = maybe False id $ do
                            { let (l,u,p,s) = lu m
                            ; let (lR,lC) = size l
                            ; let sq = abs $ lR - lC
                            ; let l' = if lC == lR then l else if lC < lR then l ||| (konst 0 (lC,sq) === ident sq) else l === (konst 0 (sq,lR) ||| ident sq)
                            ; linv <- if det l' == 0 then Nothing else Just $ inv l'
                            ; pinv <- if s == 0 then Nothing else Just $ inv p
                            ; let mat' = linv N.<> pinv N.<> m2
                            ; let matz = fromColumns $ map (toZ . roundVector) $ toColumns mat'
                            ; let cond = isUpperTri matz
                            ; let cond' = and $ map (\x -> (x - (fromIntegral $ floor x) < cutoff) || (((fromIntegral $ ceiling x) - x) < cutoff)) $ toList $ flatten mat'
                            ; return $ cond && cond'
                            }
            ; let checkAll k mns = if size (mat mns) /= size m2 then trace "Dim mismatch" False else if check $ mat mns then trace "Found one!" True else if k > ubound then trace ("Tried: "++show k) False else checkAll (k+1) (suc mns)
            ; return $ checkAll 0 n0s
            }

fromDGAMap :: StdBraid -> DGA_Map -> [Char] -> Maybe Augmentation
fromDGAMap b (DGA_Map l) chars = do
                                { l' <- mapM (\(c,a) -> (represent chars a) >>= (\vs -> return (c,vs))) l
                                ; return $ Aug b l'
                                }

compose_maps :: DGA_Map -> DGA_Map -> DGA_Map
compose_maps (DGA_Map map1) (DGA_Map map2) = DGA_Map $ (map (\(c,exp) -> (c,applyDGAMap (DGA_Map map2) exp)) map1) ++ (filter (\(c,_) -> not $ elem c $ map fst map1) map2)

applyDGAMap :: DGA_Map -> Algebra -> Algebra
applyDGAMap (DGA_Map alist) a = appmaph alist a

appmaph::[(Char,Algebra)] -> Algebra -> Algebra
appmaph [] = id
appmaph cs = plugIn (\c -> case (lookup c cs) of Just e -> e
                                                 Nothing -> G $ E c)

max_stars_repo_name: Creatorri/Legendrian-Knots-UROP

max_issues_repo_name: Creatorri/Legendrian-Knots-UROP

max_forks_repo_name: Creatorri/Legendrian-Knots-UROP