text:
{-# LANGUAGE ApplicativeDo          #-}
{-# LANGUAGE DataKinds              #-}
{-# LANGUAGE DeriveGeneric          #-}
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE GADTs                  #-}
{-# LANGUAGE KindSignatures         #-}
{-# LANGUAGE LambdaCase             #-}
{-# LANGUAGE OverloadedStrings      #-}
{-# LANGUAGE PolyKinds              #-}
{-# LANGUAGE RecordWildCards        #-}
{-# LANGUAGE ScopedTypeVariables    #-}
{-# LANGUAGE TupleSections          #-}
{-# LANGUAGE TypeApplications       #-}
{-# LANGUAGE TypeFamilies           #-}
{-# LANGUAGE TypeFamilyDependencies #-}
{-# LANGUAGE TypeOperators          #-}
{-# LANGUAGE ViewPatterns           #-}

import           Control.Arrow                         ((&&&))
import           Control.DeepSeq
import           Control.Exception
import           Control.Monad
import           Control.Monad.Trans.Class
import           Control.Monad.Trans.State.Strict
import           Control.Monad.Trans.Writer.Strict
import           Data.Bifunctor
import           Data.Either.Validation
import           Data.Finite
import           Data.Foldable
import           Data.IDX
import           Data.Kind
import           Data.List
import           Data.Monoid
import           Data.Profunctor
import           Data.Proxy
import           Data.Singletons
import           Data.Singletons.Prelude               (Sing(..))
import           Data.Singletons.TypeLits
import           Data.String
import           Data.Time.Clock
import           Data.Type.Combinator
import           Data.Type.Equality
import           Data.Type.Vector
import           GHC.Generics                          (Generic)
import           GHC.TypeLits                          as TL
import           GHC.TypeLits.Compare
import           Options.Applicative hiding            (ParserResult(..))
import           Statistics.Distribution
import           Statistics.Distribution.Uniform
import           System.Directory
import           System.FilePath
import           System.Random.MWC
import           System.Random.MWC.Distributions
import           TensorOps.Backend.BTensor
import           TensorOps.Learn.NeuralNet
import           TensorOps.Learn.NeuralNet.FeedForward
import           TensorOps.Types hiding                ((&&&))
import           Text.Printf
import           Type.Class.Higher.Util
import           Type.Family.Nat
import qualified Codec.Compression.GZip                as GZ
import qualified Control.Foldl                         as F
import qualified Data.ByteString.Lazy                  as BS
import qualified Data.Map.Strict                       as M
import qualified Data.Vector                           as V
import qualified Data.Vector.Unboxed                   as VU
import qualified Network.HTTP.Simple                   as HTTP
import qualified TensorOps.Tensor                      as TT
import qualified Text.PrettyPrint.Boxes                as B

mnistBase :: String
mnistBase = "http://yann.lecun.com/exdb/mnist"

mnistFiles :: Matrix '[N2, N2] FilePath
mnistFiles = ("train-images-idx3-ubyte" :+ "train-labels-idx1-ubyte" :+ ØV)
          :* ("t10k-images-idx3-ubyte"  :+ "t10k-labels-idx1-ubyte"  :+ ØV)
          :* ØV

data Opts = O { oRate        :: Double
              , oLayers      :: [Integer]
              , oBatch       :: Integer
              , oDataDir     :: FilePath
              , oNoConfusion :: Bool
              , oWhite       :: Bool
              , oInduce      :: Maybe (Finite 10)
              }
    deriving (Show, Eq, Generic)

opts :: Parser Opts
opts = O <$> option auto
               ( long "rate" <> short 'r' <> metavar "STEP"
              <> help "Neural network learning rate"
              <> value 0.02 <> showDefault
               )
         <*> option auto
               ( long "layers" <> short 'l' <> metavar "LIST"
              <> help "List of hidden layer sizes"
              <> value [300,100] <> showDefault
               )
         <*> option auto
               ( long "batch" <> short 'b' <> metavar "AMOUNT"
              <> help "Training batch size"
              <> value 1000 <> showDefault
               )
         <*> strOption
               ( long "data" <> short 'd' <> metavar "PATH"
              <> help "Directory to store/cache MNIST data files"
              <> value "data/mnist" <> showDefaultWith id
               )
         <*> switch
               ( long "noconfusion" <> short 'c'
              <> help "Disable confusion matrix validation and only display % error every batch"
               )
         <*> switch
               ( long "white" <> short 'w'
              <> help "Train with an eleventh \"white noise\" class to train network on negative results"
               )
         <*> optional (
               option readFin
                 ( long "induce" <> short 'i' <> metavar "DIGIT"
                <> help ("Every batch, attempt to induce an image "
                      ++ "of the given digit with the trained network"
                        )
                 )
               )
  where
    readFin :: forall n. KnownNat n => ReadM (Finite n)
    readFin = do
        i <- auto
        case packFinite i of
          Nothing -> readerError $
            printf "Number %d out of range (%d)" i (natVal (Proxy @n) - 1)
          Just x  -> return x

main :: IO ()
main = do
    O{..} <- execParser $ info (helper <*> opts)
        ( fullDesc
       <> header "tensor-ops-mnist - train neural nets on MNIST data set"
       <> progDesc (unlines ["Simple test of tensor-ops tensors (with hmatrix"
                            ,"backend) on MNIST classification challenge"
                            ]
                   )
        )

    mnistDat <- loadData oDataDir
    putStrLn "Loaded data."

    withSomeSing oWhite $ \(w :: Sing (w :: Bool)) ->
      withKnownNat (nOut w) $
        case Proxy %<=? Proxy :: 1 :<=? NOut w of
          LE Refl -> case Proxy %<=? Proxy :: 10 :<=? NOut w of
            LE Refl ->
              learn w (Proxy @(BTensorV HMatD)) mnistDat
                oRate oLayers oBatch oInduce oNoConfusion
            NLE _ -> error "impossible"
          NLE _ -> error "impossible"

loadData
    :: FilePath
    -> IO (Vec N2 [(Int, VU.Vector Int)])
loadData dataDir = do
    createDirectoryIfMissing True dataDir

    printf "Loading data from %s\n" dataDir
    mnistBs <- forM mnistFiles $ \u -> do
      BS.readFile (dataDir </> u) `catch` \(_ :: IOException) -> do
        printf "'%s' not found; downloading from %s ...\n" u mnistBase
        b <- GZ.decompress . HTTP.getResponseBody
               <$> HTTP.httpLBS (fromString (mnistBase </> u -<.> "gz"))
        BS.writeFile (dataDir </> u) b
        return b

    let f   :: Vec N2 FilePath
            -> Vec N2 BS.ByteString
            -> I (Validation [String] (FilePath, FilePath, IDXData, IDXLabels))
        f (I uIm :* I uLb :* ØV) (I bIm :* I bLb :* ØV) = I $ do
            im <- maybe (Failure [printf "Could not decode image %s." uIm]) Success
                    $ decodeIDX bIm
            lb <- maybe (Failure [printf "Could not decode labels %s." uLb]) Success
                    $ decodeIDXLabels bLb
            pure (uIm, uLb, im, lb)

        mnistDat' :: Either [String] (Vec N2 [(Int, VU.Vector Int)])
        mnistDat' = do
          imlb <- validationToEither . sequenceA $ vap f mnistFiles mnistBs
          forM imlb $ \(uIm, uLb, im, lb) ->
            maybe (Left [printf "Could not combine %s and %s." uIm uLb]) Right
              $ labeledIntData lb im

    evaluate . force
        =<< either (ioError . userError . unlines) return mnistDat'

processDat
    :: forall (n :: Nat) (l :: Nat) t.
     ( Fractional (ElemT t)
     , KnownNat n
     , KnownNat l
     , Tensor t
     )
    => (Int, VU.Vector Int)
    -> Either String (t '[n], (t '[l], Finite l))
processDat (l,d) = (,) <$> x <*> y
  where
    x :: Either String (t '[n])
    x = maybe (Left (printf "Bad input vector (Expected %d, got %d)" n (VU.length d))) Right
      . TT.fromList . map ((/255) . fromIntegral)
      $ VU.toList d
    y :: Either String (t '[l], Finite l)
    y = maybe (Left (printf "Label out of range (Got %d, expected [0,%d) )" l')) Right
      . fmap (TT.oneHot 1 0 &&& id)
      $ packFinite (fromIntegral l)
    n :: Integer
    n = natVal (Proxy @n)
    l' :: Integer
    l' = natVal (Proxy @l)

type family NOut (w :: Bool) = (n :: Nat) | n -> w where
    NOut 'False = 10
    NOut 'True  = 11

nOut
    :: Sing w
    -> Sing (NOut w)
nOut = \case
    SFalse -> SNat @10
    STrue  -> SNat @11

learn
    :: forall (t :: [Nat] -> Type) (w :: Bool) (o :: Nat).
     ( Tensor t
     , RealFloat (ElemT t)
     , NFData1 t
     , NFData (t '[784])
     , NFData (t '[o])
     , o ~ NOut w
     , 10 TL.<= o
     , 1  TL.<= o
     , KnownNat o
     )
    => Sing w
    -> Proxy t
    -> Vec N2 [(Int, VU.Vector Int)]
    -> Double
    -> [Integer]
    -> Integer
    -> Maybe (Finite 10)
    -> Bool
    -> IO ()
learn w _ dat rate layers (fromIntegral->batch) ind nc =
      withSystemRandom $ \g -> do
    dat' <- either (ioError . userError . unlines) return
          . validationToEither
          . (traverse . traverse) processDat'
          $ dat

    dat'' <- evaluate $ force dat'
    putStrLn "Data processed."

    let tXY, vXY :: [(t '[784], (t '[o], Finite o))]
        (tXY, vXY) = case dat'' of
                       I t :* I v :* ØV -> (t,v)

    net0 :: Network t 784 o
            <- genNet (layers `zip` repeat (actMap logistic)) actSoftmax g

    printf "rate: %f | batch: %d | layers: %s\n" rate batch (show layers)
    when (fromSing w) $
      putStrLn "white noise class enabled"
    forM_ ind $ \i ->
      printf "inducing: %d\n" (getFinite i)

    trainEpochs tXY ((map . second) snd vXY) g net0
  where
    processDat'
        :: (Int, VU.Vector Int)
        -> Validation [String] (t '[784], (t '[o], Finite o))
    processDat' = eitherToValidation . first (:[]) . processDat
    noiseClass :: t '[11]
    noiseClass = TT.oneHot 1 0 (finite 10)
    noiseFin   :: Finite 11
    noiseFin   = finite 10
    ind' :: Maybe (t '[o])
    ind' = TT.oneHot 1 0 . weakenN <$> ind
    trainEpochs
        :: [(t '[784], (t '[o], Finite o))]
        -> [(t '[784], Finite o)]
        -> GenIO
        -> Network t 784 o
        -> IO ()
    trainEpochs (V.fromList->tr) vd g = trainEpoch 1
      where
        trainEpoch
            :: Integer
            -> Network t 784 o
            -> IO ()
        trainEpoch e nt0 = do
            printf "[Epoch %d]\n" e
            tr' <- case w of
                SFalse -> return tr
                STrue  -> do
                  extr <- V.replicateM (V.length tr `div` 10) $ do
                            v <- genRand (uniformDistr 0 1) g
                            α <- realToFrac <$> genContVar (uniformDistr 0 1) g
                            return $ scaleT α v
                  return $ tr <> fmap (, (noiseClass, noiseFin)) extr

            queue <- evaluate . force =<< uniformShuffle tr' g

            printf "Training on %d samples in batches of %d ...\n" (V.length tr') batch

            nt1 <- trainBatch 1 queue nt0
            trainEpoch (succ e) nt1
          where
            trainBatch
                :: Integer
                -> V.Vector (t '[784], (t '[o], Finite o))
                -> Network t 784 o
                -> IO (Network t 784 o)
            trainBatch b (V.splitAt batch->(xs,xss)) nt
                | V.null xs = return nt
                | otherwise = do
              printf "Batch %d ...\n" b
              (nt', t) <- time . return $ trainAll nt ((fmap . second) fst xs)
              printf "Trained on %d samples in %s\n" (V.length xs) (show t)
              vd' :: [(t '[784], Finite o)] <- case w of
                  SFalse ->
                    return vd
                  STrue  -> do
                    extr <- replicateM (length vd `div` 10) $ do
                              v <- genRand (uniformDistr 0 1) g
                              α <- realToFrac <$> genContVar (uniformDistr 0 1) g
                              return $ scaleT α v
                    return $ vd <> fmap (, noiseFin) extr
              let tscore = F.fold (validate nt') ((fmap . second) snd xs)
                  vconf  = F.fold (confusion nt') vd'
                  confmat = ( B.vcat B.left
                                  [ B.text (printf "[%d]" (getFinite i))
                                  | i <- range :: [Finite o]
                                  ] B.<+>
                            )
                          . B.hsep 1 B.top
                          . map (\(_,r) -> B.vcat B.right
                                         . map (\(_, c) -> B.text (show c))
                                         $ M.toList r
                                )
                          $ M.toList vconf
                  vscore :: Double
                  vscore
                    | nc        = F.fold (validate nt') vd'
                    | otherwise = (/ fromIntegral (length vd')) . sum
                                . map (\(i, as) -> fromIntegral $ as M.! i)
                                $ M.toList vconf
              printf "Training:   %.2f%% error\n" ((1 - tscore) * 100)
              printf "Validation: %.2f%% error\n" ((1 - vscore) * 100)
              B.printBox confmat
              forM_ ind' $ \i -> do
                x0 <- genRand (uniformDistr 0 0.05) g
                let x1 = induceNum nt' i 1 5000 x0
                    y1 = intercalate "/"
                       . map (printf "%.2f" . realToFrac @_ @Double)
                       . TT.toList
                       $ runNetwork nt' x1
                putStrLn (renderOut x1)
                putStrLn y1

              trainBatch (succ b) xss nt'
        validate
            :: Network t 784 o
            -> F.Fold (t '[784], Finite o) Double
        validate n = (\s l -> fromIntegral s / fromIntegral l)
                 <$> lmap (uncurry f) F.sum
                 <*> F.length
          where
            f :: t '[784] -> Finite o -> Int
            f x r | TT.argMax (runNetwork n x) == r = 1
                  | otherwise                       = 0
        confusion
            :: Network t 784 o
            -> F.Fold (t '[784], Finite o) (M.Map (Finite o) (M.Map (Finite o) Integer))
        confusion n = F.Fold (\m (x, r) ->
                                let y = TT.argMax (runNetwork n x)
                                in  M.insertWith (M.unionWith (+)) y (M.singleton r 1) m
                             )
                             (M.fromList [ (i, M.fromList [ (j, 0) | j <- range ] )
                                         | i <- range
                                         ]
                               )
                             id
    trainAll
        :: Foldable f
        => Network t 784 o
        -> f (t '[784], t '[o])
        -> Network t 784 o
    trainAll = foldl' $ \nt (i,o) -> nt `deepseq`
        trainNetwork crossEntropy rate' i o nt
    rate' :: ElemT t
    rate' = realToFrac rate
    induceNum
        :: Network t 784 o
        -> t '[o]
        -> ElemT t
        -> Int
        -> t '[784]
        -> t '[784]
    induceNum n t r = go
      where
        go i x
          | i == 0    = x
          | otherwise = let x' = induceNetwork crossEntropy r t n x
                        in  x' `deepseq` go (i - 1) x'

time
    :: NFData a
    => IO a
    -> IO (a, NominalDiffTime)
time x = do
    t1 <- getCurrentTime
    y  <- evaluate . force =<< x
    t2 <- getCurrentTime
    return (y, t2 `diffUTCTime` t1)

renderOut
    :: forall t. (Tensor t, Real (ElemT t))
    => t '[784]
    -> String
renderOut = intercalate "\n"
          . ($ [])
          . appEndo
          . execWriter
          . execStateT (replicateM 28 go)
          . TT.toList
  where
    go :: StateT [ElemT t] (Writer (Endo [String])) ()
    go = do
        x <- state $ splitAt 28
        lift . tell . Endo . (++) . (:[]) $
            map (render . realToFrac)
          . concatMap (\y -> [y,y])
          $ x
    render :: Double -> Char
    render r | r <= 0.2  = ' '
             | r <= 0.4  = '.'
             | r <= 0.8  = '-'
             | r <= 1.9  = '='
             | otherwise = '#'

range
    :: forall n. KnownNat n
    => [Finite n]
range = finite <$> [0 .. natVal (Proxy @n) - 1]

max_stars_repo_name: mstksg/tensor-ops

max_issues_repo_name: mstksg/tensor-ops

max_forks_repo_name: mstksg/tensor-ops


text:
{-# LANGUAGE NoMonomorphismRestriction, Rank2Types #-}
module STMatrix (
    tests_STMatrix,
    mutatesToMatrix,
    readOnlyMatrix,

    ) where

import Control.Monad
import Control.Monad.ST
import Data.AEq
import Data.Complex( magnitude )
import Debug.Trace
import Test.Framework
import Test.Framework.Providers.QuickCheck2
import Test.QuickCheck hiding ( vector )
import qualified Test.QuickCheck as QC

import Numeric.LinearAlgebra
import qualified Numeric.LinearAlgebra.Vector as V
import qualified Numeric.LinearAlgebra.Matrix as M

import Test.QuickCheck.LinearAlgebra( TestElem(..), Dim2(..), Index2(..),
    Assocs2(..), MatrixPair(..), MatrixTriple(..) )
import qualified Test.QuickCheck.LinearAlgebra as Test

import STVector( readOnlyVector )
import Typed

tests_STMatrix = testGroup "STMatrix"
    [ testPropertyI "new_" prop_new_
    , testPropertyI "new" prop_new
    , testPropertyI "newCopy" prop_newCopy
    , testPropertyI "copyTo" prop_copyTo
    , testPropertyI "read" prop_read
    , testPropertyI "write" prop_write
    , testPropertyI "modify" prop_modify
    , testPropertyI "getElems" prop_getElems
    , testPropertyI "getElems'" prop_getElems'
    , testPropertyI "getAssocs" prop_getAssocs
    , testPropertyI "getAssocs'" prop_getAssocs'
    , testPropertyI "setElems" prop_setElems
    , testPropertyI "setAssocs" prop_setAssocs
    , testPropertyI "mapTo" prop_mapTo
    , testPropertyI "zipWithTo" prop_zipWithTo
    , testPropertyDZ "shiftDiagM_" prop_shiftDiagM_ prop_shiftDiagM_
    , testPropertyDZ "shiftDiagWithScaleToM_"
        prop_shiftDiagWithScaleM_ prop_shiftDiagWithScaleM_
    , testPropertyDZ "addTo" prop_addTo prop_addTo
    , testPropertyDZ "subTo" prop_subTo prop_subTo
    , testPropertyDZ "scaleM_" prop_scaleM_ prop_scaleM_
    , testPropertyDZ "addWithScaleM_" prop_addWithScaleM_ prop_addWithScaleM_
    , testPropertyDZ "scaleRowsM" prop_scaleRowsM_ prop_scaleRowsM_
    , testPropertyDZ "scaleColsM_" prop_scaleColsM_ prop_scaleColsM_
    , testPropertyDZ "negateTo" prop_negateTo prop_negateTo
    , testPropertyDZ "conjugateTo" prop_conjugateTo prop_conjugateTo
    , testPropertyDZ "swapRows" prop_swapRows prop_swapRows
    , testPropertyDZ "swapCols" prop_swapCols prop_swapCols
    ]

prop_new_ t (Dim2 n) =
    (M.dim $ typed t $ M.create $ M.new_ n) === n

prop_new t (Dim2 n) e =
    (M.create $ M.new n e) === (typed t $ M.constant n e)

prop_newCopy t x =
    (M.create $ M.newCopy x) === x
  where
    _ = typed t x

prop_copyTo t (MatrixPair x y) = runST $
    x `readOnlyMatrix` \mx ->
    y `mutatesToMatrix` x $ \my ->
        M.copyTo my mx
  where
    _ = typed t x

prop_read t (Index2 n i) =
    forAll (typed t `fmap` Test.matrix n) $ \x -> runST $
        x `readOnlyMatrix` \mx -> do
            e <- M.read mx i
            return $ e === M.at x i

prop_write t (Index2 n i) e =
    forAll (Test.matrix n) $ \x -> runST $
        x `mutatesToMatrix` (x `M.update` [(i,e)]) $ \mx -> do
            M.write mx i e
  where
    _ = e == t

prop_modify t (Index2 n i) (Blind f) =
    forAll (Test.matrix n) $ \x -> runST $
        x `mutatesToMatrix`
            (typed t $ x `M.update` [(i, f $ M.at x i)]) $ \mx ->
                M.modify mx i f

prop_getElems t x =
    forAll arbitrary $ \(Blind f) -> runST $
        x `mutatesToMatrix` (M.map f x) $ \mx -> do
            es <- M.getElems mx
            M.mapTo mx f mx
            return $ es === M.elems (M.map f x)
  where
    _ = typed t x

prop_getElems' t x =
    forAll arbitrary $ \(Blind f) -> runST $
        x `mutatesToMatrix` (M.map f x) $ \mx -> do
            es <- M.getElems' mx
            M.mapTo mx f mx
            return $ es === M.elems x
  where
    _ = typed t x

prop_getAssocs t x =
    forAll arbitrary $ \(Blind f) -> runST $
        x `mutatesToMatrix` (M.map f x) $ \mx -> do
            ies <- M.getAssocs mx
            M.mapTo mx f mx
            return $ ies === M.assocs (M.map f x)
  where
    _ = typed t x

prop_getAssocs' t x =
    forAll arbitrary $ \(Blind f) -> runST $
        x `mutatesToMatrix` (M.map f x) $ \mx -> do
            ies <- M.getAssocs' mx
            M.mapTo mx f mx
            return $ ies === M.assocs x
  where
    _ = typed t x

prop_setElems t x =
    forAll (QC.vector $ m*n) $ \es -> runST $
        x `mutatesToMatrix` (M.fromList (m,n) es) $ \mx ->
            M.setElems mx es
  where
    (m,n) = M.dim x
    _ = typed t x

prop_setAssocs t (Assocs2 mn ies) =
    forAll (Test.matrix mn) $ \x -> runST $
        x `mutatesToMatrix` (typed t $ M.update x ies) $ \mx ->
            M.setAssocs mx ies

prop_mapTo t (Blind f) = binaryProp t
    (\x -> M.map f x)
    (\dst mx -> M.mapTo dst f mx)

prop_zipWithTo t (Blind f) = ternaryProp t
    (\x y -> M.zipWith f x y)
    (\dst mx my -> M.zipWithTo dst f mx my)

prop_shiftDiagM_ t a =
    forAll (Test.vector (min m n)) $ \s -> runST $
        s `readOnlyVector` \ms ->
        a `mutatesToMatrix` (M.shiftDiag s a) $ \ma ->
            M.shiftDiagM_ ms ma
  where
    (m,n) = M.dim a
    _ = typed t a

prop_shiftDiagWithScaleM_ t e a =
    forAll (Test.vector (min m n)) $ \s -> runST $
        s `readOnlyVector` \ms ->
        a `mutatesToMatrix` (M.shiftDiagWithScale e s a) $ \ma ->
            M.shiftDiagWithScaleM_ e ms ma
  where
    (m,n) = M.dim a
    _ = typed t a

prop_addTo t = ternaryProp t M.add M.addTo

prop_subTo t = ternaryProp t M.sub M.subTo

prop_scaleM_ t e a = runST $
    a `mutatesToMatrix` (M.scale e a) $ \ma ->
        M.scaleM_ e ma
  where
    _ = typed t a

prop_addWithScaleM_ t alpha (MatrixPair x y) = runST $
    x `readOnlyMatrix` \mx ->
    y `mutatesToMatrix` (M.addWithScale alpha x y) $ \my ->
        M.addWithScaleM_ alpha mx my
  where
    _ = typed t x

prop_scaleRowsM_ t a =
    forAll (Test.vector m) $ \s -> runST $
        s `readOnlyVector` \ms ->
        a `mutatesToMatrix` (M.scaleRows s a) $ \ma ->
            M.scaleRowsM_ ms ma
  where
    (m,n) = M.dim a
    _ = typed t a

prop_scaleColsM_ t a =
    forAll (Test.vector n) $ \s -> runST $
        s `readOnlyVector` \ms ->
        a `mutatesToMatrix` (M.scaleCols s a) $ \ma ->
            M.scaleColsM_ ms ma
  where
    (m,n) = M.dim a
    _ = typed t a

prop_negateTo t = binaryProp t M.negate M.negateTo

prop_conjugateTo t = binaryProp t M.conjugate M.conjugateTo

prop_swapRows t =
    forAll arbitrary $ \(Index2 (m,n) (i1,_)) ->
    forAll (choose (0,m-1)) $ \i2 ->
    forAll (typed t `fmap` Test.matrix (m,n)) $ \a ->
        let a' = M.update a $ [ ((i2,j),e) | (j,e) <- V.assocs (M.row a i1)
                              ] ++
                              [ ((i1,j),e) | (j,e) <- V.assocs (M.row a i2)
                              ]
        in runST $
            a `mutatesToMatrix` a' $ \ma ->
            M.swapRows ma i1 i2

prop_swapCols t =
    forAll arbitrary $ \(Index2 (m,n) (_,j1)) ->
    forAll (choose (0,n-1)) $ \j2 ->
    forAll (typed t `fmap` Test.matrix (m,n)) $ \a ->
        let a' = M.update a $ [ ((i,j2),e) | (i,e) <- V.assocs (M.col a j1)
                              ] ++
                              [ ((i,j1),e) | (i,e) <- V.assocs (M.col a j2)
                              ]
        in runST $
            a `mutatesToMatrix` a' $ \ma ->
            M.swapCols ma j1 j2

binaryProp :: (AEq e, Arbitrary e, Show e, Storable e, Testable a)
            => e
            -> (Matrix e -> Matrix e)
            -> (forall s . STMatrix s e -> STMatrix s e -> ST s a)
            -> Property
binaryProp t imm_f f = let
    prop1 t x = runST $
        x `mutatesToMatrix` (imm_f x) $ \mx -> f mx mx
      where _ = typed t x
    prop2 t (MatrixPair x y) = runST $ let _ = typed t x in
        x `readOnlyMatrix` \mx ->
        y `mutatesToMatrix` (imm_f x) $ \my -> f my mx
      where _ = typed t x
    in (    label "1" prop1
        .&. label "2" prop2
       )

ternaryProp :: (AEq e, Arbitrary e, Show e, Storable e, Testable a)
            => e
            -> (Matrix e -> Matrix e -> Matrix e)
            -> (forall s . STMatrix s e -> STMatrix s e -> STMatrix s e -> ST s a)
            -> Property
ternaryProp t imm_f f = let
    prop1 t x = runST $
        x `mutatesToMatrix` (imm_f x x) $ \mx -> f mx mx mx
      where _ = typed t x
    prop2a t (MatrixPair x y) = runST $
        x `readOnlyMatrix` \mx ->
        y `mutatesToMatrix` (imm_f x y) $ \my -> f my mx my
      where _ = typed t x
    prop2b t (MatrixPair x y) = runST $
        x `mutatesToMatrix` (imm_f x y) $ \mx ->
        y `readOnlyMatrix` \my -> f mx mx my
      where _ = typed t x
    prop3 t (MatrixTriple x y z) = runST $
        x `readOnlyMatrix` \mx ->
        y `readOnlyMatrix` \my ->
        z `mutatesToMatrix` (imm_f x y) $ \mz -> f mz mx my
      where _ = typed t x
    in (    label "1" prop1
        .&. label "2a" prop2a
        .&. label "2b" prop2b
        .&. label "3" prop3
       )

readOnlyMatrix :: (Storable e, AEq e, Testable prop, Show e)
               => Matrix e
               -> (STMatrix s e ->  ST s prop)
               -> ST s Property
readOnlyMatrix x = mutatesToMatrix x x

mutatesToMatrix :: (Storable e, AEq e, Testable prop, Show e)
                => Matrix e
                -> Matrix e
                -> (STMatrix s e -> ST s prop)
                -> ST s Property
mutatesToMatrix x x_new f = do
    mx <- M.newCopy x
    prop <- f mx
    x' <- M.freeze mx
    return $ prop .&. (x' === x_new)

ext: hs

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
module Statistics.ANOVA (anova, fisherLSD) where
import Data.List
import Statistics.Distribution hiding (mean, variance)
import Statistics.Distribution.FDistribution
import Statistics.Distribution.StudentT
import Statistics.Sample
import qualified Data.Vector as V

{-# INLINE treatments #-}
treatments :: (Ord a,Eq a) => V.Vector (a, Double) -> V.Vector a
treatments input = V.fromList $! nub $! V.toList $! V.map fst input --not good

treatmentMeans :: (Ord a, Eq a) => V.Vector a -> V.Vector (a, Double) -> V.Vector (a,Double)
treatmentMeans treatments input = V.map (\x -> x input) meanFuncs
  where meanFuncs = V.map treatmentMean $! treatments

treatmentMean :: (Ord a, Eq a) => a -> V.Vector (a,Double) -> (a,Double)
treatmentMean treatment input = (treatment, avg)
  where avg = mean $! V.map snd $! filterTreatment treatment $! input

treatmentSize treatment input = V.length $! filterTreatment treatment $! input

{-# INLINE filterTreatment #-}
filterTreatment :: (Ord a, Eq a) => a -> V.Vector (a, Double) -> V.Vector (a, Double)
filterTreatment treatment = V.filter (\(x,y) -> x == treatment)

sampleMean input = mean $! V.map snd input

ssTreatment :: (Ord a, Eq a) => V.Vector a -> V.Vector (a, Double) -> Double
ssTreatment treatments input = V.sum $! V.zipWith (\size mean -> size*(mean-sMean)*(mean-sMean)) tSizes tMeans
  where sMean = sampleMean input
        tMeans = V.map snd $! treatmentMeans treatments input
        tSizes = V.map (\x -> fromIntegral $! x input) $! V.map treatmentSize $! treatments --input

ssError :: (Ord a, Eq a) => V.Vector (a, Double) -> Double
ssError input = V.sum $! V.map (\x -> x input) $! V.map treatmentError $! treatments input

treatmentError treatment input = V.sum $! V.map (\x -> (x-tMean)*(x-tMean)) filteredInput
  where filteredInput = V.map snd $! filterTreatment treatment input
        tMean = snd $! treatmentMean treatment input

fValue :: (Ord a, Eq a) => V.Vector a -> V.Vector (a, Double) -> Double
fValue treatments input = ((ssTreatment treatments input)/(ssError input)) * ((sampleSize - numberOfTreatments)/(numberOfTreatments -1))
  where numberOfTreatments = fromIntegral $! V.length $! treatments
        sampleSize = fromIntegral $! V.length input

fDist treatments input = fDistribution (numberOfTreatments-1) (sampleSize-numberOfTreatments)
  where numberOfTreatments = fromIntegral $! V.length treatments
        sampleSize = fromIntegral $! V.length input

pValue input = 1 - cumulative (fDist t input) (fValue t input)
  where t = treatments input

anova pvalue input
  | pValue input < pvalue = True
  | otherwise = False

tTest :: (Ord a, Eq a) => Double -> a -> a -> V.Vector (a, Double) -> Bool
tTest pvalue treatment1 treatment2 inputs = abs (tTestT treatment1 treatment2 inputs) > critical
  where critical = quantile (studentT (degreesOfFreedom treatment1 treatment2 inputs)) (1 - pvalue/2)

tTestT :: (Ord a, Eq a) => a -> a -> V.Vector (a, Double) -> Double
tTestT treatment1 treatment2 inputs = (mean1-mean2)/sqrt(var1/size1+var2/size2)
  where mean1 = snd $! treatmentMean treatment1 inputs
        mean2 = snd $! treatmentMean treatment2 inputs
        size1 = fromIntegral $! V.length $! filterTreatment treatment1 inputs
        size2 = fromIntegral $! V.length $! filterTreatment treatment2 inputs
        var1 = variance $! V.map snd $! filterTreatment treatment1 inputs
        var2 = variance $! V.map snd $! filterTreatment treatment2 inputs

degreesOfFreedom
  :: (Ord a, Eq a) => a -> a -> V.Vector (a, Double) -> Double
degreesOfFreedom treatment1 treatment2 inputs = numerator / denominator
  where size1 = fromIntegral $! V.length $! filterTreatment treatment1 inputs
        size2 = fromIntegral $! V.length $! filterTreatment treatment2 inputs
        var1 = variance $! V.map snd $! filterTreatment treatment1 inputs
        var2 = variance $! V.map snd $! filterTreatment treatment2 inputs
        numerator = (var1/size1+var2/size2)^2
        denominator = (((var1/size1)^2)/(size1-1)) + (((var2/size2)^2)/(size2-1))

noDifferenceInMean pvalue treatment input = (:) treatment $!
                                            map fst $!
                                            filter (\(x,y) -> not y) $!
                                            map (\x -> (x, tTest pvalue treatment x input)) othertreatments
  where othertreatments = delete treatment $! V.toList $! treatments input

equivalenceList pvalue input = nub $! V.toList $! V.map (\x -> sort $! noDifferenceInMean pvalue x input) $! treatments input

fisherLSD :: (Eq a, Ord a) => Double -> V.Vector (a, Double) -> [[a]]
fisherLSD pvalue input = case (anova pvalue input) of
  False -> [V.toList $! treatments input]
  True -> equivalenceList pvalue input

max_stars_repo_name: richardfergie/statistics-anova

max_issues_repo_name: richardfergie/statistics-anova

max_forks_repo_name: richardfergie/statistics-anova


text:
{-# LANGUAGE FlexibleContexts, FlexibleInstances, UndecidableInstances #-}

-- Pure Haskell implementation of standard BFGS algorithm with cubic
-- line search, with the BFGS step exposed to provide access to
-- intermediate Hessian estimates.

-- Author: Ian Ross, for OpenBrain Ltd

module Fuml.Optimisation.BFGS
       ( BFGSOpts (..)
       , LineSearchOpts (..)
       , bfgs, bfgsWith
       , bfgsInit
       , bfgsStep, bfgsStepWith )
where

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Devel
import Data.Default
--import Debug.Trace

-- Type synonyms mostly just for documentation purposes...
type Point = Vector Double
type Gradient = Vector Double
type Direction = Vector Double
type Fn = Point -> Double
type GradFn = Point -> Vector Double
type Hessian = Matrix Double

--------------------------------------------------------------------------------
--
--  NaN CHECKING
--

class Nanable a where
  hasnan :: a -> Bool

instance Nanable Double where
  hasnan = isNaN

instance (Nanable a, Container c a) => Nanable (c a) where
  hasnan = not . null . (find hasnan)

--------------------------------------------------------------------------------
--
--  BFGS
--

-- Options for BFGS solver: point tolerance, gradient tolerance,
-- maximum iterations.
--
data BFGSOpts = BFGSOpts { ptol :: Double
                         , gtol :: Double
                         , maxiters :: Int } deriving Show

instance Default BFGSOpts where
  def = BFGSOpts 1.0E-7 1.0E-7 200

-- BFGS solver data to carry between steps: current point, function
-- value at point, gradient at point, current direction, current
-- Hessian estimate, maximum line search step.
--
data BFGS = BFGS { p :: Point
                 , fp :: Double
                 , g :: Gradient
                 , xi :: Direction
                 , h :: Matrix Double
                 , stpmax :: Double } deriving Show

-- Main solver interface with default options.
--
bfgs :: Fn -> GradFn -> Point -> Either String (Point, Hessian)
bfgs f df p0 = bfgsWith def f df p0 (ident $ size p0)

{- Collection solver state into an infinite list: useful as "take n $
-- bfgsCollect f df b0"
--
bfgsCollect :: Fn -> GradFn -> BFGS -> [BFGS]
bfgsCollect f df b0 = map snd $ iterate step (False,b0)
  where step (cvg,b) = if cvg then b else bfgsStep f df b

-}
-- Utility function to set up initial BFGS state for bfgsCollect.
--
bfgsInit :: Fn -> GradFn -> Point -> Either String BFGS
bfgsInit f df p0 = case (hasnan f0, hasnan g0) of
  (False, False) -> Right $ BFGS p0 f0 g0 (-g0) (ident n) (maxStep p0)
  _ -> Left $ nanMsg p0 (Just f0) (Just g0)
  where n = size p0 ; f0 = f p0 ; g0 = df p0

-- Main iteration routine: sets up initial BFGS state, then steps
-- until converged or maximum iterations exceeded.
--
bfgsWith :: BFGSOpts -> Fn -> GradFn -> Point -> Hessian -> Either String (Point, Hessian)
bfgsWith opt@(BFGSOpts _ _ maxiters') f df p0 h0 =
  case (hasnan f0, hasnan g0) of
    (False, False) -> go 0 b0
    errs -> Left $ nanMsg p0 (Just f0) (Just g0)
  where go iters b =
          if iters > maxiters'
          then Left "maximum iterations exceeded in bfgs"
          else case bfgsStepWith opt f df b of
            Left err -> Left err
            Right (True, b') -> Right (p b', h b')
            Right (False, b') -> go (iters+1) $ {- trace (show $ (iters, fp b')) -} b'
        f0 = f p0 ; g0 = df p0
        b0 = BFGS p0 f0 g0 (-g0) h0 (maxStep p0)

-- Do a BFGS step with default parameters.
--
bfgsStep :: Fn -> GradFn -> BFGS -> Either String (Bool, BFGS)
bfgsStep = bfgsStepWith def

-- Main BFGS step routine.  This is a more or less verbatim
-- translation of the description in Section 10.7 of Numerical Recipes
-- in C, 2nd ed.
--
bfgsStepWith :: BFGSOpts -> Fn -> GradFn -> BFGS -> Either String (Bool, BFGS)
bfgsStepWith (BFGSOpts ptol' gtol' _) f df (BFGS p fp g xi h stpmax) =
  case lineSearch f p fp g xi stpmax of
    Left err -> Left err
    Right (pn, fpn) ->
      if hasnan gn
      then Left $ nanMsg pn Nothing (Just gn)
      else if cvg
           then Right (True, BFGS pn fpn gn xi h stpmax)
           else Right (False, BFGS pn fpn gn xin hn stpmax)
      where gn = df pn ; dp = pn - p ; dg = gn - g ; hdg = h #> dg
            dpdg = dp `dot` dg ; dghdg = dg `dot` hdg
            hn = h + ((dpdg + dghdg) / dpdg^2) `scale` (dp `outer` dp) -
                 (1/dpdg) `scale` (h <> (dg `outer` dp) + (dp `outer` dg) <> h)
            xin = -hn #> gn
            cvg = maxabsratio dp p < ptol' || maxabsratio' (fpn `max` 1) gn p < gtol'

-- Generate error messages for NaN production in function and gradient
-- calculations.
--
nanMsg :: Point -> Maybe Double -> Maybe Gradient -> String
nanMsg p' fval grad = "NaNs produced: p = " ++ show p' ++
                     maybe "" (("  fval = " ++) . show) fval ++
                     maybe "" (("  grad = " ++) . show) grad

--------------------------------------------------------------------------------
--
--  LINE SEARCH
--

-- Options for line search routine: point tolerance, "acceptable
-- decrease" parameter.
--
data LineSearchOpts = LineSearchOpts { xtol :: Double
                                     , alpha :: Double } deriving Show

instance Default LineSearchOpts where
  def = LineSearchOpts 1.0E-7 1.0E-4

-- Maximum line search step length.
--
maxStep :: Point -> Double
maxStep p0 = (100 * (norm_2 p0 `max` n))
  where n = fromIntegral (size p0)

-- Line search with default parameters.
--
lineSearch :: Fn -> Point -> Double -> Gradient -> Direction ->
              Double -> Either String (Point, Double)
lineSearch = lineSearchWith def

-- Main line search routine.  This is kind of nasty to translate into
-- functional form because of the switching about between the
-- quadratic and cubic approximations.  It works, but it could be
-- prettier.
--
lineSearchWith :: LineSearchOpts -> Fn -> Point -> Double -> Gradient ->
                  Direction -> Double -> Either String (Point, Double)
lineSearchWith (LineSearchOpts xtol alpha) func xold fold g pin stpmax =
  go 1.0 Nothing
  where p = if pinnorm > stpmax then (stpmax/pinnorm) `scale` pin else pin
        pinnorm = norm_2 pin
        slope = g `dot` p
        lammin = xtol / (maxabsratio p xold)

        go :: Double -> Maybe (Double,Double) -> Either String (Point,Double)
        go lam pass = if hasnan fnew
                      then Left $ nanMsg xnew (Just fnew) Nothing
                      else case check xnew fnew of
                        Just xandf -> Right xandf
                        Nothing ->
                          case pass of
                            -- First time.
                            Nothing -> go (lambound $ quadlam fnew) $ Just (lam,fnew)
                            -- Subsequent times.
                            Just val2 -> case cubiclam fnew val2 of
                              Right newlam -> go (lambound newlam) $ Just (lam,fnew)
                              Left err -> Left err
          where xnew = xold + lam `scale` p
                fnew = func xnew

                -- Check for convergence or a "sufficiently large" step.
                check :: Vector Double -> Double -> Maybe (Vector Double,Double)
                check x f =
                  if lam < lammin then Just (xold,fold)
                  else if f <= fold + alpha * lam * slope
                       then Just (x,f)
                       else Nothing

                -- Keep step length within bounds.
                lambound lam' = max (0.1 * lam) (min lam' (0.5 * lam))

                -- Quadratic and cubic approximations to better step
                -- value.
                quadlam fnew = -slope / (2 * (fnew - fold - slope))
                cubiclam fnew (lam2,f2) =
                  if a == 0
                  then Right (-slope / (2 * b))
                  else if disc < 0
                       then Left "Roundoff problem in lineSearch"
                       else Right $ (-b + sqrt disc) / (3 * a)
                    where rhs1 = fnew - fold - lam * slope
                          rhs2 = f2 - fold - lam2 * slope
                          a = (rhs1 / lam^2 - rhs2 / lam2^2) / (lam - lam2)
                          b = (-lam2 * rhs1 / lam^2 + lam * rhs2 / lam2^2) /
                              (lam - lam2)
                          disc = b^2 - 3 * a * slope

-- Utility functions for ratio testing.
--
absratio :: Double -> Double -> Double
absratio n d = abs n / (abs d `max` 1)

absratio' :: Double -> Double -> Double -> Double
absratio' scale n d = abs n / (abs d `max` 1) / scale

maxabsratio :: Vector Double -> Vector Double -> Double
maxabsratio n d = maxElement $ zipVectorWith absratio n d

maxabsratio' :: Double -> Vector Double -> Vector Double -> Double
maxabsratio' scale n d = maxElement $ zipVectorWith (absratio' scale) n d

--------------------------------------------------------------------------------
--
--  TEST FUNCTIONS
--

-- Simple test function.  Minimum at (3,4):
--
-- *BFGS> bfgs tstf1 tstgrad1 (fromList [-10,-10])
-- Right (fromList [3.0,4.0])
--

{-
tstf1 :: Fn
tstf1 p = (x-3)^2 + (y-4)^2
  where [x,y] = toList p

tstgrad1 :: GradFn
tstgrad1 p = fromList [2*(x-3),2*(y-4)]
  where [x,y] = toList p

-- Rosenbrock's function.  Minimum at (1,1):
--
-- *BFGS> bfgs tstf2 tstgrad2 (fromList [-10,-10])
-- Right (fromList [0.9999999992103538,0.9999999985219549])

tstf2 :: Fn
tstf2 p = (1-x)^2 + 100*(y-x^2)^2
  where [x,y] = toList p

tstgrad2 :: GradFn
tstgrad2 p = fromList [2*(x-1) - 400*x*(y-x^2), 200*(y-x^2)]
  where [x,y] = toList p

-- Test function from Numerical Recipes.  Minimum at (-2.0,0.89442719):
--
-- *BFGS> bfgs tstfnr tstgradnr nrp0
-- Right (fromList [-1.9999999564447526,0.8944271925873616])

tstfnr :: Fn
tstfnr p = 10*(y^2*(3-x)-x^2*(3+x))^2+(2+x)^2/(1+(2+x)^2)
  where [x,y] = toList p

tstgradnr :: GradFn
tstgradnr p = fromList [20*(y^2*x3m-x^2*x3p)*(-y^2-6*x-3*x^2)+
                        2*x2p/(1+x2p^2)-2*x2p^3/(1+x2p^2)^2,
                        40*(y^2*x3m-x^2*x3p)*y*x3m]
  where [x,y] = toList p
        x3m = 3 - x
        x3p = 3 + x
        x2p = 2 + x

nrp0 :: Point
nrp0 = fromList [0.1,4.2]

-- Test function to check NaN handling.
--
-- *BFGS> bfgs nantstf nantstgrad (fromList [-10,-10])
-- Left "function application returned NaN"

nantstf :: Fn
nantstf p = log x + (x-3)^2 + (y-4)^2
  where [x,y] = toList p

nantstgrad :: GradFn
nantstgrad p = fromList [1/x+2*(x-3),2*(y-4)]
  where [x,y] = toList p
-}

max_stars_repo_name: avctrh/filopodia

max_issues_repo_name: avctrh/filopodia

max_forks_repo_name: avctrh/filopodia


text:
{-# LANGUAGE BangPatterns, DeriveDataTypeable, DeriveGeneric,
    FlexibleContexts #-}
-- |
-- Module    : Statistics.Sample.Powers
-- Copyright : (c) 2009, 2010 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Very fast statistics over simple powers of a sample.  These can all
-- be computed efficiently in just a single pass over a sample, with
-- that pass subject to stream fusion.
--
-- The tradeoff is that some of these functions are less numerically
-- robust than their counterparts in the 'Statistics.Sample' module.
-- Where this is the case, the alternatives are noted.

module Statistics.Sample.Powers
    (
    -- * Types
      Powers

    -- * Constructor
    , powers

    -- * Descriptive functions
    , order
    , count
    , sum

    -- * Statistics of location
    , mean

    -- * Statistics of dispersion
    , variance
    , stdDev
    , varianceUnbiased

    -- * Functions over central moments
    , centralMoment
    , skewness
    , kurtosis

    -- * References
    -- $references
    ) where

import Control.Monad.ST
import Data.Aeson            (FromJSON, ToJSON)
import Data.Binary           (Binary(..))
import Data.Data             (Data, Typeable)
import Data.Vector.Binary    ()
import Data.Vector.Unboxed   ((!))
import GHC.Generics          (Generic)
import Numeric.SpecFunctions (choose)
import Prelude hiding (sum)
import Statistics.Function   (indexed)
import qualified Data.Vector          as V
import qualified Data.Vector.Generic  as G
import qualified Data.Vector.Storable as SV
import qualified Data.Vector.Unboxed  as U
import qualified Data.Vector.Unboxed.Mutable as MU
import qualified Statistics.Sample.Internal  as S

newtype Powers = Powers (U.Vector Double)
    deriving (Eq, Read, Show, Typeable, Data, Generic)

instance FromJSON Powers
instance ToJSON Powers

instance Binary Powers where
    put (Powers v) = put v
    get = fmap Powers get

-- | O(/n/) Collect the /n/ simple powers of a sample.
--
-- Functions computed over a sample's simple powers require at least a
-- certain number (or /order/) of powers to be collected.
--
-- * To compute the /k/th 'centralMoment', at least /k/ simple powers
--   must be collected.
--
-- * For the 'variance', at least 2 simple powers are needed.
--
-- * For 'skewness', we need at least 3 simple powers.
--
-- * For 'kurtosis', at least 4 simple powers are required.
--
-- This function is subject to stream fusion.
powers :: G.Vector v Double =>
          Int                   -- ^ /n/, the number of powers, where /n/ >= 2.
       -> v Double
       -> Powers
powers k sample
  | k < 2     = error "Statistics.Sample.powers: too few powers"
  | otherwise = runST $ do
      acc <- MU.replicate l 0
      G.forM_ sample $ \x ->
        let loop !i !xk
              | i == l    = return ()
              | otherwise = do MU.write acc i . (+ xk) =<< MU.read acc i
                               loop (i+1) (xk * x)
        in loop 0 1
      fmap Powers $ U.unsafeFreeze acc
  where
    l = k + 1
{-# SPECIALIZE powers :: Int -> U.Vector  Double -> Powers #-}
{-# SPECIALIZE powers :: Int -> V.Vector  Double -> Powers #-}
{-# SPECIALIZE powers :: Int -> SV.Vector Double -> Powers #-}

-- | The order (number) of simple powers collected from a 'sample'.
order :: Powers -> Int
order (Powers pa) = U.length pa - 1

-- | Compute the /k/th central moment of a sample.  The central
-- moment is also known as the moment about the mean.
centralMoment :: Int -> Powers -> Double
centralMoment k p@(Powers pa)
    | k < 0 || k > order p =
                  error ("Statistics.Sample.Powers.centralMoment: "
                         ++ "invalid argument")
    | k == 0    = 1
    | otherwise = (/n) . S.sum . U.map go . indexed . U.take (k+1) $ pa
  where
    go (i , e) = (k `choose` i) * ((-m) ^ (k-i)) * e
    n = U.head pa
    m = mean p

-- | Maximum likelihood estimate of a sample's variance.  Also known
-- as the population variance, where the denominator is /n/.  This is
-- the second central moment of the sample.
--
-- This is less numerically robust than the variance function in the
-- 'Statistics.Sample' module, but the number is essentially free to
-- compute if you have already collected a sample's simple powers.
--
-- Requires 'Powers' with 'order' at least 2.
variance :: Powers -> Double
variance = centralMoment 2

-- | Standard deviation.  This is simply the square root of the
-- maximum likelihood estimate of the variance.
stdDev :: Powers -> Double
stdDev = sqrt . variance

-- | Unbiased estimate of a sample's variance.  Also known as the
-- sample variance, where the denominator is /n/-1.
--
-- Requires 'Powers' with 'order' at least 2.
varianceUnbiased :: Powers -> Double
varianceUnbiased p@(Powers pa)
    | n > 1     = variance p * n / (n-1)
    | otherwise = 0
  where n = U.head pa

-- | Compute the skewness of a sample. This is a measure of the
-- asymmetry of its distribution.
--
-- A sample with negative skew is said to be /left-skewed/.  Most of
-- its mass is on the right of the distribution, with the tail on the
-- left.
--
-- > skewness . powers 3 $ U.to [1,100,101,102,103]
-- > ==> -1.497681449918257
--
-- A sample with positive skew is said to be /right-skewed/.
--
-- > skewness . powers 3 $ U.to [1,2,3,4,100]
-- > ==> 1.4975367033335198
--
-- A sample's skewness is not defined if its 'variance' is zero.
--
-- Requires 'Powers' with 'order' at least 3.
skewness :: Powers -> Double
skewness p = centralMoment 3 p * variance p ** (-1.5)

-- | Compute the excess kurtosis of a sample.  This is a measure of
-- the \"peakedness\" of its distribution.  A high kurtosis indicates
-- that the sample's variance is due more to infrequent severe
-- deviations than to frequent modest deviations.
--
-- A sample's excess kurtosis is not defined if its 'variance' is
-- zero.
--
-- Requires 'Powers' with 'order' at least 4.
kurtosis :: Powers -> Double
kurtosis p = centralMoment 4 p / (v * v) - 3
    where v = variance p

-- | The number of elements in the original 'Sample'.  This is the
-- sample's zeroth simple power.
count :: Powers -> Int
count (Powers pa) = floor $ U.head pa

-- | The sum of elements in the original 'Sample'.  This is the
-- sample's first simple power.
sum :: Powers -> Double
sum (Powers pa) = pa ! 1

-- | The arithmetic mean of elements in the original 'Sample'.
--
-- This is less numerically robust than the mean function in the
-- 'Statistics.Sample' module, but the number is essentially free to
-- compute if you have already collected a sample's simple powers.
mean :: Powers -> Double
mean p@(Powers pa)
    | n == 0    = 0
    | otherwise = sum p / n
    where n     = U.head pa

-- $references
--
-- * Besset, D.H. (2000) Elements of statistics.
--   /Object-oriented implementation of numerical methods/
--   ch. 9, pp. 311–331.
--   <http://www.elsevier.com/wps/product/cws_home/677916>
--
-- * Anderson, G. (2009) Compute /k/th central moments in one
--   pass. /quantblog/. <http://quantblog.wordpress.com/2009/02/07/compute-kth-central-moments-in-one-pass/>


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE FlexibleContexts #-}
module STC.CompletionFieldR2S1 where

import           Array.UnboxedArray          as AU
import           Control.Monad               as M
import           Control.Monad.Parallel      as MP
import           Control.Parallel.Strategies
import           Data.Array.Repa             as R
import           Data.Complex
import           Data.List                   as L
import           Data.Vector.Storable        as VS
import           Data.Vector.Unboxed         as VU
import           DFT.Plan
import           Filter.Utils
import           Image.IO
import           Image.Transform
import           System.FilePath
import           System.Random
import           Text.Printf
import           Types
import           Utils.Array
import           Utils.Time
import           FokkerPlanck.FourierSeries
import FokkerPlanck.Analytic
import qualified Numeric.LinearAlgebra as NL
import qualified Numeric.LinearAlgebra.HMatrix as NL
import            STC.Utils
import FokkerPlanck.BrownianMotion (thetaPlus)
import Utils.Diagram

-- Implementation of [Williams and Jacbos, 1997]

{-# INLINE checkN #-}
checkN :: Int -> Int -> Int
checkN maxN n
  | n < 0 = checkN maxN (n + maxN)
  | n >= maxN = checkN maxN (n - maxN)
  | otherwise = n

{-# INLINE rotateST #-}
rotateST ::
     (R.Source r Double)
  => R.Array r DIM3 Double
  -> Int
  -> R.Array U DIM3 Double
rotateST arr 0 = computeUnboxedS . delay $ arr
rotateST arr n' =
  let (Z :. nf :. nx :. ny) = extent arr
      n = checkN nf n'
      deg = (-2) * pi / (fromIntegral nf) * fromIntegral n
  in rotate25D deg (fromIntegral $ center nx, fromIntegral $ center ny) $
     R.backpermute
       (extent arr)
       (\(Z :. k :. i :. j) -> (Z :. (mod (k + nf - n) nf) :. i :. j))
       arr

{-# INLINE timeReversal #-}
timeReversal :: (R.Source r e) => Array r DIM3 e -> Array D DIM3 e
timeReversal arr =
  let (Z :. nf :. _ :. _) = extent arr
      n = div nf 2
  in R.backpermute
       (extent arr)
       (\(Z :. k :. i :. j) -> (Z :. (mod (k + n) nf) :. i :. j))
       arr

{-# INLINE rotateSTTR #-}
rotateSTTR ::
     (R.Source r Double)
  => R.Array r DIM3 Double
  -> R.Array U DIM3 Double
rotateSTTR arr =
  let (Z :. _ :. nx :. ny) = extent arr
  in rotate25D pi (fromIntegral $ center nx, fromIntegral $ center ny) arr

{-# INLINE makeR2S1Plan #-}
makeR2S1Plan :: (R.Source r e) => DFTPlan -> R.Array r DIM3 e -> IO DFTPlan
makeR2S1Plan oldPlan arr = do
  let (Z :. orientations :. rows :. cols) = extent arr
  vec3D <-
    VS.map (:+ 0) . VS.fromList <$>
    M.replicateM (orientations * rows * cols) randomIO
  lock <- getFFTWLock
  fst <$>
    (dft1dGPlan lock oldPlan [orientations, rows, cols] [1, 2] vec3D >>= \(plan, vec) ->
       idft1dGPlan lock plan [orientations, rows, cols] [1, 2] vec)

{-# INLINE convolve #-}
convolve ::
     (R.Source r (Complex Double))
  => DFTPlan
  -> Array U DIM3 Double
  -> Array r DIM2 (Complex Double)
  -> IO (VU.Vector Double)
convolve dftPlan arr3D arr2DF = do
  let (Z :. orientations :. rows :. cols) = extent arr3D
      planID3D = DFTPlanID DFT1DG [orientations, rows, cols] [1, 2]
      inversePlanID = DFTPlanID IDFT1DG [orientations, rows, cols] [1, 2]
  vec3DF <-
    dftExecute dftPlan planID3D . VU.convert . VU.map (:+ 0) . toUnboxed $ arr3D
  let arr3DF =
        fromUnboxed (Z :. orientations :. rows :. cols) . VS.convert $ vec3DF
  vec <-
    dftExecute
      dftPlan
      inversePlanID
      (VU.convert . toUnboxed . computeS . R.traverse2 arr3DF arr2DF const $ \f3 f2 idx@(Z :. k :. i :. j) ->
         f3 idx * f2 (Z :. i :. j))
  return . VU.map realPart . VS.convert $ vec

{-# INLINE shareWeightST #-}
shareWeightST ::
     (R.Source r Double)
  => DFTPlan
  -> Array U DIM3 Double
  -> Array r DIM3 Double
  -> IO (Array D DIM3 Double)
shareWeightST dftPlan arr arrG = do
  let (Z :. nf :. _ :. _) = extent arrG
      (Z :. orientations :. rows :. cols) = extent arr
  arrF <-
    fmap (fromUnboxed (extent arr) . VS.convert) .
    dftExecute dftPlan (DFTPlanID DFT1DG [orientations, rows, cols] [1, 2]) .
    VU.convert . VU.map (:+ 0) . toUnboxed $
    arr
  xs <-
    MP.mapM
      (\i ->
         convolve
           dftPlan
           (computeS . makeFilter2D . rotateST arrG $ i)
           (R.slice arrF $ (Z :. i :. All :. All)))
      [0 .. nf - 1]
  return . delay . fromUnboxed (extent arrG) . L.foldl1' (VU.zipWith (+)) $ xs

{-# INLINE computeInitialDistribution #-}
computeInitialDistribution ::
     Int -> Int -> Int -> Double -> [R2S1RPPoint] -> R.Array U DIM3 Double
computeInitialDistribution xLen yLen numOrientation delta xs =
  let deltaTheat = 2 * pi / fromIntegral numOrientation
      xShift = div xLen 2
      (xMin, xMax) =
        if odd xLen
          then (-xShift, xShift)
          else (-xShift, xShift - 1)
      yShift = div yLen 2
      (yMin, yMax) =
        if odd yLen
          then (-yShift, yShift)
          else (-yShift, yShift - 1)
      vec =
        toUnboxedVector .
        AU.accum (+) 0 ((0, xMin, yMin), (numOrientation - 1, xMax, yMax)) .
        L.map
          (\(R2S1RPPoint (x, y, theta, _)) ->
             ( ((floor $ theta * pi / 180 / deltaTheat), round (fromIntegral x / delta), round (fromIntegral y / delta))
             , (1 / (fromIntegral . L.length $ xs)))) $
        xs
   in fromUnboxed (Z :. numOrientation :. xLen :. yLen) vec

-- Implementation of [Williams and Zweck, 2003]

{-# INLINE computeBias #-}
computeBias ::
     Int -> Int -> Int -> [R2S1RPPoint] -> R.Array D DIM3 Double
computeBias xLen yLen numOrientation xs =
  let xShift = div xLen 2
      (xMin, xMax) =
        if odd xLen
          then (-xShift, xShift)
          else (-xShift, xShift - 1)
      yShift = div yLen 2
      (yMin, yMax) =
        if odd yLen
          then (-yShift, yShift)
          else (-yShift, yShift - 1)
      vec =
        toUnboxedVector .
        AU.accum (+) 0 ((xMin, yMin), (xMax, yMax)) .
        L.map (\(R2S1RPPoint (x, y, _, _)) -> ((x, y), 1)) $
        xs
   in R.traverse
        (fromUnboxed (Z :. xLen :. yLen) vec)
        (const (Z :. numOrientation :. xLen :. yLen)) $ \f (Z :. _ :. i :. j) ->
        f (Z :. i :. j)

{-# INLINE computeInitialEigenVec #-}
computeInitialEigenVec ::
     Int -> Int -> Int -> [R2S1RPPoint] -> R.Array D DIM3 Double
computeInitialEigenVec xLen yLen numOrientation xs =
  let xShift = div xLen 2
      (xMin, xMax) =
        if odd xLen
          then (-xShift, xShift)
          else (-xShift, xShift - 1)
      yShift = div yLen 2
      (yMin, yMax) =
        if odd yLen
          then (-yShift, yShift)
          else (-yShift, yShift - 1)
      vec =
        toUnboxedVector .
        AU.accum (+) 0 ((xMin, yMin), (xMax, yMax)) .
        L.map
          (\(R2S1RPPoint (x, y, _, _)) ->
             ((x, y), (1 / (fromIntegral . L.length $ xs)))) $
        xs
   in R.traverse
        (fromUnboxed (Z :. xLen :. yLen) vec)
        (const (Z :. numOrientation :. xLen :. yLen)) $ \f (Z :. _ :. i :. j) ->
        f (Z :. i :. j) / fromIntegral numOrientation

{-# INLINE computeInitialEigenVecGaussian #-}
computeInitialEigenVecGaussian ::
     Int -> Int -> Int -> [(Int, Int, Double)] -> R.Array D DIM3 Double
computeInitialEigenVecGaussian xLen yLen numOrientation xs =
  let xShift = div xLen 2
      (xMin, xMax) =
        if odd xLen
          then (-xShift, xShift)
          else (-xShift, xShift - 1)
      yShift = div yLen 2
      (yMin, yMax) =
        if odd yLen
          then (-yShift, yShift)
          else (-yShift, yShift - 1)
      vec =
        toUnboxedVector .
        AU.accum (+) 0 ((xMin, yMin), (xMax, yMax)) .
        L.map (\(x, y, v) -> ((x, y), v)) $
        xs
  in R.traverse
       (fromUnboxed (Z :. xLen :. yLen) vec)
       (const (Z :. numOrientation :. xLen :. yLen)) $ \f (Z :. _ :. i :. j) ->
       f (Z :. i :. j)

{-# INLINE eigenVectorR2S1 #-}
eigenVectorR2S1 ::
     (R.Source s Double, R.Source s1 Double, R.Source s2 Double)
  => DFTPlan
  -> FilePath
  -> R.Array s DIM3 Double
  -> R.Array s DIM3 Double
  -> Int
  -> Bool
  -> R.Array s1 DIM3 Double
  -> R.Array s2 DIM3 Double
  -> IO (R.Array D DIM3 Double)
eigenVectorR2S1 plan folderPath filterSource filterSink n writeFlag bias inputSource = do
  let sourceDist = R.zipWith (*) bias inputSource
      -- s = R.sumAllS sourceDist
      s = L.maximum . R.toList $ sourceDist
      sourceVec = computeS $ R.map (/ s) sourceDist
      (Z :. numOrientation :. _ :. _) = extent inputSource
  printCurrentTime (show n)
  sourceEigenVec <- shareWeightST plan sourceVec filterSource
  -- plotThetaDimension folderPath (printf "R2S1_%d_" n) (31,-4) .
  --   R.backpermute
  --     (extent sourceEigenVec)
  --     (\(Z :. k :. i :. j) ->
  --        (Z :.
  --         (mod (k + numOrientation - (div numOrientation 2)) numOrientation :: Int) :.
  --         i :.
  --         j)) $
  --   sourceEigenVec
  -- sinkEigenVec <-
  --   timeReversal <$>
  --   shareWeightST
  --     plan
  --     (computeS . R.zipWith (*) bias $ sourceEigenVec')
  --     filterSink
  -- let sourceEigenVec =
  --       R.zipWith (\x y -> x + 0.2 * y) sourceEigenVec' sinkEigenVec
  when
    writeFlag
    (plotImageRepa
       (folderPath </> printf "Source_%d.png" n)
       (ImageRepa 8 .
        computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.map (\(x) -> (abs x) ** (1/3)) . R.sumS . rotate3D $
        sourceEigenVec))
  return sourceEigenVec

powerMethod ::
     (R.Source s1 Double, R.Source s2 Double)
  => DFTPlan
  -> FilePath
  -> (R.Array U DIM3 Double)
  -> Int
  -> Bool
  -> String
  -> Double
  -> (R.Array s1 DIM3 Double)
  -> (R.Array s2 DIM3 Double)
  -> IO (R.Array U DIM3 Double)
powerMethod oldPlan folderPath filterSource numIteration writeFlag idStr threshold bias' initialEigenVec = do
  let (Z :. numOrientation :. xLen :. yLen) = extent filterSource
      filterSink = rotateSTTR filterSource
      bias = delay bias'
  plotImageRepa
    (folderPath </> "Bias.png")
    (ImageRepa 8 .
     computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.sumS . rotate3D $
     bias)
  plan <- makeR2S1Plan oldPlan filterSource
  source <-
    M.foldM
      (\input n ->
         eigenVectorR2S1
           plan
           folderPath
           filterSource
           filterSink
           n
           writeFlag
           bias
           input)
      (delay initialEigenVec)
      [1 .. numIteration]
  -- print . R.toList . R.slice source $ (Z :. All :. (div xLen 2) :. (div yLen 2))
  -- plotThetaDimension folderPath "R2S1_Source_" (31,-4) .
  --   R.backpermute
  --     (extent source)
  --     (\(Z :. k :. i :. j) ->
  --        (Z :.
  --         (mod (k + numOrientation - (div numOrientation 2)) numOrientation :: Int) :.
  --         i :.
  --         j)) $
  --   source
  sink <- shareWeightST plan (rotateSTTR $ bias *^ source) filterSink
  plotImageRepa
    (folderPath </> "Sink.png")
    (ImageRepa 8 .
     computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.sumS . rotate3D $
     sink)
  -- plotThetaDimension folderPath "R2S1_Sink_" (31,-4) .
  --   R.backpermute
  --     (extent sink)
  --     (\(Z :. k :. i :. j) ->
  --        (Z :.
  --         (mod (k + numOrientation - (div numOrientation 2)) numOrientation :: Int) :.
  --         i :.
  --         j)) $
  --   sink
  let completion = computeS $ R.zipWith (*) source sink
      completionR2 = R.sumS . rotate3D $ completion
      completionR2Vec = toUnboxed completionR2
      completionMax = VU.maximum completionR2Vec
      completionMin = VU.minimum completionR2Vec
  plotImageRepa
    (folderPath </> printf "Completion%s.png" idStr)
    (ImageRepa 8 .
     computeS .
     R.extend (Z :. (1 :: Int) :. All :. All) .
     R.map (\x -> sqrt $ (x - completionMin) / (completionMax - completionMin)) $
     completionR2)
  -- plotThetaDimension folderPath "R2S1_Completion_" (31, -4) .
  --   R.backpermute
  --     (extent completion)
  --     (\(Z :. k :. i :. j) ->
  --        (Z :.
  --         (mod (k + numOrientation - (div numOrientation 2)) numOrientation :: Int) :.
  --         i :.
  --         j)) $
  --   completion
  return completion

eigenVectorR2S1Analytic ::
     Int -> Double -> Double -> Double -> Double -> [(Double, Double)] -> IO [[Double]]
eigenVectorR2S1Analytic !oris !sigma !tau !gamma !delta !xs = do
  let !numPoints = L.length xs
      !locArr = fromListUnboxed (Z :. L.length xs) xs
      !deltaTheta = 2 * pi / fromIntegral oris
      -- !weights =
      --   [1 / 16, 1 / 8, 1 / 16, 1 / 8, 1 / 4, 1 / 8, 1 / 16, 1 / 8, 1 / 16]
      -- !origins =
      --   L.zipWith
      --     (\w (i, j) -> (i * delta, j * delta, w))
      --     weights
      --     [(i, j) | i <- [-1 .. 1], j <- [-1 .. 1]]
      transitionMatrixArr' =
        R.traverse locArr (\_ -> (Z :. numPoints :. oris :. numPoints :. oris)) $ \f (Z :. i1 :. j1 :. i2 :. j2) ->
          let rowIdx = i1 * oris + j1
              colIdx = i2 * oris + j2
          in if rowIdx /= colIdx && i1 /= i2
             -- then L.foldl'
             --        (\s (i, j, w) ->
             --           s +
             --           w *
             --           computePji
             --             sigma
             --             tau
             --             (R2S1RP
             --                (fst . f $ (Z :. i2))
             --                (snd . f $ (Z :. i2))
             --                (deltaTheta * fromIntegral j2)
             --                gamma)
             --             (R2S1RP
             --                (i + (fst . f $ (Z :. i1)))
             --                (j + (snd . f $ (Z :. i1)))
             --                (deltaTheta * fromIntegral j1)
             --                gamma))
             --        0
             --        origins
             then computePji
                    sigma
                    tau
                    (R2S1RP
                       ((fst . f $ (Z :. i2)))
                       ((snd . f $ (Z :. i2)))
                       (deltaTheta * fromIntegral j2)
                       gamma)
                    (R2S1RP
                       (fst . f $ (Z :. i1))
                       (snd . f $ (Z :. i1))
                       (deltaTheta * fromIntegral j1)
                       gamma)
             else 0
  transitionMatrixArr <- computeUnboxedP transitionMatrixArr'
      -- transitionMatrixArr = computeUnboxedS transitionMatrixArr'
  let !size = numPoints * oris
      transitionMatrix =
        deepSeqArray transitionMatrixArr . (size NL.>< size) . R.toList $
        transitionMatrixArr
      (eigVal, eigVec) = NL.eig transitionMatrix
      pairs =
        L.reverse .
        L.sortOn (realPart . fst) .
        L.filter
          (\(c, _) ->
             let (m, p) = polar c
             in abs (p) < 1e-10) $
        L.zip (NL.toList eigVal) (L.map NL.toList . NL.toColumns $ eigVec)
      eigenVector =
        fromListUnboxed (Z :. numPoints :. oris) .
        L.map magnitude . snd . L.head $
        pairs
  print . L.map (polar . fst) . L.take 10 $ pairs
  print . L.map (polar . fst) . L.take 10 . L.reverse $ pairs
  print . L.maximum . L.map realPart . snd . L.head $ pairs
  print . L.minimum . L.map realPart . snd . L.head $ pairs
  return . L.map (\i -> R.toList . R.slice eigenVector $ (Z :. i :. All)) $
    [0 .. numPoints - 1]

computeContourR2S1Analytic ::
     (R.Source s Double)
  => FilePath
  -> (R.Array U DIM3 Double)
  -> Int
  -> Int
  -> Int
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> R.Array s DIM3 Double
  -> [(Double, Double)]
  -> IO (R.Array U DIM3 Double)
computeContourR2S1Analytic !folderPath !filterSource !oris !rows !cols !sigma !tau !gamma !std !delta !bias !xs' = do
  let xs =
        L.map
          (\(x', y') -> (fromIntegral . round $ x', fromIntegral . round $ y'))
          xs'
  eigenVector <-
    eigenVectorR2S1Analytic
      oris
      sigma
      tau
      gamma
      delta
      (L.map (\(a, b) -> (a * delta, b * delta)) xs)
  let ys =
        L.concat $
        L.zipWith
          (\(x, y) thetas ->
             L.zipWith
               (\k theta -> ((k, round x, round y), theta))
               [0 ..]
               thetas)
          xs
          eigenVector
      (!minR, !maxR) = computeRange rows
      (!minC, !maxC) = computeRange cols
      eigenSourceArr =
        fromUnboxed (Z :. oris :. cols :. rows) .
        toUnboxedVector .
        AU.accum (+) 0 ((0, minC, minR), (oris - 1, maxC, maxR)) $
        ys
      eigenSinkArr = rotateSTTR eigenSourceArr
      completionArr = computeUnboxedS $ eigenSourceArr *^ eigenSinkArr
      filterSink = rotateSTTR filterSource
  plotR2S1Array
    (folderPath </> "EigenSource.eps")
    (fromIntegral rows)
    (fromIntegral cols)
    (fromIntegral rows / 16)
    xs
    eigenSourceArr
  plotR2S1Array
    (folderPath </> "EigenSink.eps")
    (fromIntegral rows)
    (fromIntegral cols)
    (fromIntegral rows / 16)
    xs
    eigenSinkArr
  plotR2S1Array
    (folderPath </> "EigenCompletion.eps")
    (fromIntegral rows)
    (fromIntegral cols)
    (fromIntegral rows / 16)
    xs
    completionArr
  plan <- makeR2S1Plan emptyPlan filterSource
  source <- shareWeightST plan eigenSourceArr filterSource
  plotImageRepa
    (folderPath </> "Source.png")
    (ImageRepa 8 .
     computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.sumS . rotate3D $
     source)
  plotR2S1Array
    (folderPath </> "Source.eps")
    (fromIntegral rows)
    (fromIntegral cols)
    (fromIntegral rows / 16)
    xs
    source
  sink <- shareWeightST plan eigenSinkArr filterSink
  plotImageRepa
    (folderPath </> "Sink.png")
    (ImageRepa 8 .
     computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.sumS . rotate3D $
     sink)
  plotR2S1Array
    (folderPath </> "Sink.eps")
    (fromIntegral rows)
    (fromIntegral cols)
    (fromIntegral rows / 16)
    xs
    sink
  let completion = computeS $ R.zipWith (*) source sink
      completionR2 = R.sumS . rotate3D $ completion
      completionR2Vec = toUnboxed completionR2
      completionMax = VU.maximum completionR2Vec
      completionMin = VU.minimum completionR2Vec
  print (completionMin, completionMax)
  plotImageRepa
    (folderPath </> printf "Completion.png")
    (ImageRepa 8 .
     computeS .
     R.extend (Z :. (1 :: Int) :. All :. All) .
     -- R.map (\x -> (abs x) ** (1 / 2)) .
     R.map (\x -> sqrt $ (x - completionMin) / (completionMax - completionMin)) $
     completionR2)
  plotR2S1Array
    (folderPath </> "Completion.eps")
    (fromIntegral rows)
    (fromIntegral cols)
    (fromIntegral rows / 16)
    xs
    completion
  return completion

computeContourR2S1Tangent ::
     FilePath
  -> (R.Array U DIM3 Double)
  -> Int
  -> Int
  -> Int
  -> Double
  -> [(Double, Double)]
  -> IO (R.Array U DIM3 Double)
computeContourR2S1Tangent !folderPath !filterSource !oris !rows !cols !delta !xs = do
  let (cMin, cMax) = computeRange cols
      (rMin, rMax) = computeRange rows
      !deltaTheta = 2 * pi / fromIntegral oris
      initArr =
        fromUnboxed (Z :. oris :. cols :. rows) .
        toUnboxedVector .
        AU.accum (+) 0 ((0, cMin, rMin), (oris - 1, cMax, rMax)) .
        L.concatMap
          (\(x, y) ->
             let phi = atan2 y x
                 idx1 = floor $ ((phi `thetaPlus` (pi / 2)) + pi) / deltaTheta
                 idx2 = floor $ ((phi `thetaPlus` ((-pi) / 2)) + pi) / deltaTheta
             in [((idx1, round x, round y), 1), ((idx2, round x, round y), 1)]) $
        xs
      filterSink = rotateSTTR filterSource
  plan <- makeR2S1Plan emptyPlan filterSource
  source <- shareWeightST plan initArr filterSource
  plotImageRepa
    (folderPath </> "Source.png")
    (ImageRepa 8 .
     computeS .
     R.extend (Z :. (1 :: Int) :. All :. All) .
     R.map (\x -> (abs x) ** (1 / 3)) . R.sumS . rotate3D $
     source)
  sink <- shareWeightST plan initArr filterSink
  plotImageRepa
    (folderPath </> "Sink.png")
    (ImageRepa 8 .
     computeS .
     R.extend (Z :. (1 :: Int) :. All :. All) .
     R.map (\x -> (abs x) ** (1 / 3)) . R.sumS . rotate3D $
     sink)
  let completion = computeS $ R.zipWith (*) source sink
      completionR2 = R.sumS . rotate3D $ completion
      completionR2Vec = toUnboxed completionR2
      completionMax = VU.maximum completionR2Vec
      completionMin = VU.minimum completionR2Vec
  print (completionMin, completionMax)
  plotImageRepa
    (folderPath </> printf "Completion.png")
    (ImageRepa 8 .
     computeS .
     R.extend (Z :. (1 :: Int) :. All :. All) .
     R.map (\x -> (abs x) ** (1 / 2)) .
     R.map (\x -> sqrt $ (x - completionMin) / (completionMax - completionMin)) $
     completionR2)
  return completion

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE FlexibleContexts #-}

module Taylor where

import Data.Complex

-- also this is a vector space
data Taylor a = Integ a (Taylor a)

instance Show a => Show (Taylor a) where
  show = show . take 10 . toList

toList (Integ x xs) = x:toList xs

zero = Integ 0 zero

instance Num a => Num (Taylor a) where
  (a `Integ` as) + (b `Integ` bs)  =  (a + b)  `Integ`  (as + bs)
  (a `Integ` as) * (b `Integ` bs) = (a*b) `Integ`  (as * (b `Integ` bs) + (a `Integ` as) * bs)
  fromInteger x = integ (fromInteger x) 0
  negate (a `Integ` as) = negate a `Integ` negate as

constant :: Num a => a -> Taylor a
constant x = integ x 0

integ :: a -> Taylor a -> Taylor a
integ = Integ

sinx :: Taylor Integer
sinx = integ 0 cosx

cosx :: Taylor Integer
cosx = integ 1 (-sinx)

expx :: Taylor (Complex Double)
expx = integ 1 expx

-- expkx :: Taylor Double -> Taylor Double
-- expkx k = integ 1 (k*expx)

varx :: Taylor (Complex Double)
varx = integ 0 1

(*^) :: Num t => t -> Taylor t -> Taylor t
a *^ (Integ k ks) = Integ (a * k) (a *^ ks)

evalPart :: Fractional p => [p] -> p -> p -> p
evalPart [] _ _  = 0
evalPart (k:ks) n y = k + evalPart (((y/n) *) <$> ks) (n+1) y

eval :: Fractional p => Taylor p -> Int -> p -> p
eval t n x = evalPart (take n $ toList t) 1 x

-- f . g
-- integ ((f . g) 0) (D (f . g))
-- integ (f (g 0)) (D f . g * g')

compose :: (Eq a, Num a) => Taylor a -> Taylor a -> Taylor a
compose (Integ f0 f') g@(Integ g0 g') = Integ fg0 ((f' `compose` g) * g')
  where
    -- fg0 is the evaluation of f at point g0. (An infinite sum in the general case, but...)
    fg0 = if g0 == 0 then f0 else error "fg0: computation diverges"

-- >>> varx
-- [0,1,0,0,0,0,0,0,0,0]

-- >>> 3 * varx
-- [0.0,3.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]

-- >>> expx `compose` ((constant (0 :+ 1) * varx))
-- [1.0 :+ 0.0,
--  0.0 :+ 1.0,
--  (-1.0) :+ 0.0,
--  0.0 :+ (-1.0),
--  1.0 :+ 0.0,
--  0.0 :+ 1.0,
--  (-1.0) :+ 0.0,
--  0.0 :+ (-1.0),
--  1.0 :+ 0.0,
--  0.0 :+ 1.0]

max_stars_repo_name: nicolabotta/DSLsofMath

max_issues_repo_name: nicolabotta/DSLsofMath

max_forks_repo_name: nicolabotta/DSLsofMath


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

{-# LANGUAGE AllowAmbiguousTypes   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedLabels      #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-|
Module      : Grenade.Layers.Transpose
Description : Layer that transposes its input
Maintainer  : Theo Charalambous
License     : BSD2
Stability   : experimental
-}

module Grenade.Layers.Transpose
  (
  -- * Layer Definition
    Transpose (..)

  -- * Helper function
  , initTranspose
  ) where

import           Control.DeepSeq                   (NFData (..))

import           Data.Kind                         (Type)
import           Data.Maybe                        (fromJust)
import           Data.Proxy
import           Data.Serialize
import           GHC.TypeLits

import qualified Numeric.LinearAlgebra             as LA
import           Numeric.LinearAlgebra.Static      (R)
import qualified Numeric.LinearAlgebra.Static      as H

import           Grenade.Core
import           Grenade.Layers.Internal.Transpose
import           Grenade.Onnx

 -- todo: we can probably use a type family to represent this much better
data Transpose :: Nat
               -> Shape
               -> Shape
               -> Type where
  Transpose  :: ( KnownNat dimensions )
          => !(R dimensions)
          -> Transpose dimensions input output

instance Show (Transpose dimensions input output) where
  show (Transpose mat) = "Transpose " ++ show mat

instance UpdateLayer (Transpose dimensions input output) where
  type Gradient (Transpose dimensions input output) = ()
  runUpdate _ x _  = x
  reduceGradient _ = ()

instance ( KnownNat dimensions ) => RandomLayer (Transpose dimensions input output) where
  createRandomWith _ _ = pure initTranspose

-- | Initialize a tranpose layer that is equivalent to the identity function
initTranspose :: forall dimensions input output. ( KnownNat dimensions )
              => Transpose dimensions input output
initTranspose =
  let ds    = fromIntegral $ natVal (Proxy :: Proxy dimensions)
      perms = H.fromList [1..ds] :: R dimensions
  in Transpose perms

instance ( KnownNat dimensions ) => Serialize (Transpose dimensions input output) where
  put (Transpose perms) = putListOf put . LA.toList . H.extract $ perms
  get                   = do
    perms <- maybe (fail "Vector of incorrect size") return . H.create . LA.fromList =<< getListOf get
    return $ Transpose perms

instance ( KnownNat i, KnownNat j, KnownNat k, KnownNat l, KnownNat a, KnownNat b, KnownNat c, KnownNat d )
  => Layer (Transpose 4 ('D4 i j k l) ('D4 a b c d)) ('D4 i j k l) ('D4 a b c d) where

  type Tape (Transpose 4 ('D4 i j k l) ('D4 a b c d)) ('D4 i j k l) ('D4 a b c d) = ()

  runForwards (Transpose perms) (S4D x)
    = let n  = fromIntegral $ natVal (Proxy :: Proxy i)
          c  = fromIntegral $ natVal (Proxy :: Proxy j)
          h  = fromIntegral $ natVal (Proxy :: Proxy k)
          w  = fromIntegral $ natVal (Proxy :: Proxy l)
          x' = H.extract x
          perms' = H.extract perms
          r  = transpose4d [n, c, h, w] perms' x'
      in  ((), S4D . fromJust . H.create $ r)

  runBackwards = undefined

instance OnnxOperator (Transpose dimensions input output) where
  onnxOpTypeNames _ = ["Transpose"]

instance OnnxLoadable (Transpose 4 input output) where
  loadOnnxNode _ node = readIntsAttribute "perm" node >>= formatPerm
    where
      formatPerm :: [Int] -> Either OnnxLoadFailure (Transpose 4 input output)
      formatPerm [_, _, n, c, h, w]
        = let centered = map (\x -> fromIntegral $ x - 2) [n, c, h, w]
              perms :: R 4 = H.fromList centered
          in  return $ Transpose perms
      formatPerm _ = loadFailureReason "Permutation shape incorrect of Transpose"

instance NFData (Transpose dims input output) where
  rnf (Transpose perms) = rnf perms

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
-- |
-- Module      : BattleHack.Render
-- Description :
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental
-- Portability : POSIX

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC pragmas
--------------------------------------------------------------------------------------------------------------------------------------------

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module BattleHack.Render where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Control.Monad (forM_, liftM, when)
import Control.Lens
import Control.Applicative
import Data.Complex
import Data.Fixed (mod')
import qualified Data.Set as S

import qualified Graphics.Rendering.Cairo as Cairo

import BattleHack.Types
import BattleHack.Lenses
import BattleHack.Utilities.Vector
import BattleHack.Utilities.General
import BattleHack.Utilities.Math
import qualified BattleHack.Piano as Piano

--------------------------------------------------------------------------------------------------------------------------------------------
-- Data
--------------------------------------------------------------------------------------------------------------------------------------------

--------------------------------------------------------------------------------------------------------------------------------------------
-- Functions
--------------------------------------------------------------------------------------------------------------------------------------------
-- General rendering functions -------------------------------------------------------------------------------------------------------------
-- |
polygon :: [Complex Double] -> Cairo.Render ()
polygon (p:oints) = vectorise Cairo.moveTo p >> forM_ (oints) (vectorise Cairo.lineTo) >> Cairo.closePath

-- Piano -----------------------------------------------------------------------------------------------------------------------------------
-- |
-- TODO: Create utility functions for significant key layout points
key :: PianoSettings -> Int -> Cairo.Render ()
key piano i = do
  polygon . map (+Piano.keyorigin piano i) $ keylayout
  let (r, g, b, a) = colour i in Cairo.setSourceRGBA r g b a
  Cairo.fill

  keylabel piano i
  Cairo.fill
  where
    keylayout = Piano.layout (piano-->keysize) (piano-->indent) (piano-->mid) (Piano.keylayout i) :: [Complex Double]
    colour ikey
      | (piano-->keys) !! ikey            = (0.3, 0.12, 0.22, 1.0)
      | Just ikey == (piano-->active)     = (0.3, 0.12, 0.22, 1.0)
      | mod ikey 12 `elem` Piano.naturals = (0.4, 1/7 * fromIntegral (ikey `mod` 7), 0.75, 1.00)
      | otherwise                         = (0.0, 0.0,                               0.00, 1.00)

-- |
keylabel :: PianoSettings -> Int -> Cairo.Render ()
keylabel piano i = do
  Cairo.setFontSize (if (i `mod` 12) `elem` Piano.naturals then 48 else 22)
  Cairo.selectFontFace "Old English" Cairo.FontSlantNormal Cairo.FontWeightNormal
  Cairo.setSourceRGBA 0.20 0.12 0.08 1.00
  centredText (Piano.keyorigin piano i + o + dotwise (*) sz (0.5:+0.90)) (Piano.notenameFromKeyIndex i)
  where
    (o, sz) = Piano.keybounds piano (Piano.keylayout i)

-- |
claviature :: PianoSettings -> Cairo.Render ()
claviature settings = do
  forM_ (zipWith const [0..] (settings-->keys)) $ \i -> do
    key settings i

--------------------------------------------------------------------------------------------------------------------------------------------
-- |
sinewave :: AppState -> Cairo.Render ()
sinewave appstate = do
  vectorise Cairo.moveTo o
  forM_ [o + (θ:+(48*sin (fromIntegral n/ fromIntegral fps')*sin (θ*0.04))) | θ <- [0,5..600]] (vectorise Cairo.lineTo)
  Cairo.setLineWidth 12
  Cairo.setSourceRGBA 0.62 0.94 0.44 1.0
  Cairo.stroke
  where
    o    = 6.32 * 50:+50
    n    = appstate-->animation.frame
    fps' = appstate-->animation.fps

--------------------------------------------------------------------------------------------------------------------------------------------
-- |
debugHUD :: AppState -> Cairo.Render ()
debugHUD appstate = do
  Cairo.setFontSize 14
  Cairo.selectFontFace "Helvetica" Cairo.FontSlantNormal Cairo.FontWeightNormal
  Cairo.setSourceRGBA 0.05 0.59 0.62 1.00

  let keywidth = appstate-->piano.keysize.real
  when (Piano.insideLeft (appstate-->piano) $ (mouse' & real %~ flip mod' keywidth)) $ do
    centredText (100:+30) "Inside left"

  when (Piano.insideRight (appstate-->piano) $ (mouse' & real %~ flip mod' keywidth)) $ do
    centredText (100:+45) "Inside right"

  maybe pass (\i -> centredText (100:+60) $ Piano.notenameFromKeyIndex i) (appstate-->piano.active)
  -- when (Piano.insideLeft (appstate-->piano) $ (mouse' & real %~ flip mod' keywidth)) $ do
    -- centredText (appstate-->piano.origin + (100:+30)) "Inside left"

  where
    mouse' = appstate-->inputstate.mouse

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

-- |
textsize :: String -> Cairo.Render Vector
textsize text = do
  extents <- Cairo.textExtents text
  return $ Cairo.textExtentsWidth extents :+ Cairo.textExtentsHeight extents

-- | Borrowed from Southpaw
-- TODO: General anchor (?)
centredText :: Vector -> String -> Cairo.Render ()
centredText centre text = do
  size <- textsize text
  vectorise Cairo.moveTo $ centre - 0.5*size
  Cairo.showText text

-- Experiments  ----------------------------------------------------------------------------------------------------------------------------

-- |
radial :: AppState -> Cairo.Render ()
radial appstate = Cairo.withRadialPattern a b c d e f $ \ pattern -> do
  Cairo.patternAddColorStopRGBA pattern 0 1 1 1 1
  Cairo.patternAddColorStopRGBA pattern 1 0 0 0 1
  Cairo.setSource pattern
  Cairo.arc mx my 76.8 0 (2*π)
  Cairo.fill

  Cairo.arc mx my 76.8 0 (2*π)
  Cairo.setSourceRGBA 0.0 0.0 0.0 1.0
  Cairo.setLineWidth  12
  Cairo.stroke
  where
    (mx:+my) = appstate-->inputstate.mouse
    [a, b, c, d, e, f] = [115.2, 102.4, 25.6, 102.4, 102.4, 128.0]

-- Menu ------------------------------------------------------------------------------------------------------------------------------------

-- |
overlay :: AppState -> Cairo.Render ()
overlay appstate = do
  -- TODO: Writing animation
  -- Cairo.setOperator Cairo.OperatorDestIn
  Cairo.setFontSize 180
  Cairo.selectFontFace ("Verdana" :: String) Cairo.FontSlantNormal Cairo.FontWeightBold
  Cairo.setSourceRGBA 1.0 1.0 1.0 0.8
  -- centredText ((winx:+winy) * 0.5) "CHORDIAL"
  -- Cairo.setLineWidth 12
  size' <- textsize "CUERDIAL"
  vectorise Cairo.moveTo $ (winsize - size') * 0.5
  Cairo.textPath "CUERDIAL"
  Cairo.clip
  -- Cairo.stroke
  where
    winsize = dotwise (*) (appstate-->piano.keysize) (7:+1) + 2*(appstate-->piano.origin) -- TODO: Find canvas size directly

-- Menu ------------------------------------------------------------------------------------------------------------------------------------


text:
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE FlexibleContexts #-}
module Statistics.BBVI.Gradient
  ( gradientScore
  , gradientReparam
  , DistInvariant(..)
  , Samples
  )
where

import qualified Data.Vector                   as V
import           Statistics.BBVI.Propagator     ( DistCell(..)
                                                , Gradient
                                                , Memory
                                                , dist
                                                )
import           Statistics.BBVI.Class

type Samples = V.Vector

-- | invariant data for a variational distribution cell
data DistInvariant a = DistInvariant { weight :: !Double
                                     -- ^ total occurences of z in log
                                     -- likelihood
                                     , prior :: !a
                                     -- ^ prior on latent variable
                                     , rhoF :: !(Gradient ->
                                                 DistCell a ->
                                                 (Memory, V.Vector Double))
                                     -- ^ calculates change in memory,
                                     -- and stepsize
                                     }

-- | internal helper function for building gradient propagators
gradient
  :: (DistUtil a)
  => (DistCell a -> Double -> c -> Double -> V.Vector Double)
  -> DistInvariant a
  -> DistCell a
  -> (Double, V.Vector Double, Samples c)
  -> DistCell a3
gradient f (DistInvariant {..}) no (nFactors, like, samples) = U
  memory'
  (V.zipWith (*) rho' gr)
 where
  summed =
    V.foldl' (V.zipWith (+)) (V.replicate (nParams (dist no)) 0.0)
      $ V.zipWith (f no nFactors) samples like
  gr              = fmap (/ (fromIntegral $ V.length samples)) summed
  (memory', rho') = rhoF gr no

-- | helper function for transforming log joint to score gradient
-- propagator
gradientScore
  :: Dist a c
  => DistInvariant a
  -> DistCell a -- ^ current variational distribution
  -> (Double, V.Vector Double, Samples c)   -- ^ n factors in log
  -- joint calculation, log joint calculation on each sample of \(z\),
  -- samples
  -> DistCell a -- ^ update to variational distribution
gradientScore gp@(DistInvariant {..}) = gradient f gp
 where
  f (Node _time _memory d) nFactors s l = fmap
    (* (l + nFactors / weight * (logProb prior s - logProb d s)))
    (paramGradOfLogQ d s)
  f (U _ _) _ _ _ = error "called gradient update on update cell"

-- | helper function for transforming log joint to reparameterization
-- gradient propagator
gradientReparam
  :: (Differentiable a Double)
  => DistInvariant a
  -> DistCell a -- ^ current variational distribution
  -> (Double, V.Vector Double, Samples Double)
  -- ^ n factors in log joint calculation, gradient of joint wrt \(z\)
  -- for each sample, un-transformed samples
  -> DistCell a -- ^ update to variational distribution
gradientReparam gp@(DistInvariant {..}) = gradient f gp
 where
  f (Node _time _memory d) nFactors s g = fmap
    (* ( g
       + nFactors
       / weight
       * ( sampleGradOfLogQ prior (transform d s)
         - sampleGradOfLogQ d     (transform d s)
         )
       )
    )
    (gradTransform d s)
  f (U _ _) _ _ _ = err

err :: a
err = error "called gradient update on update cell"

max_stars_repo_name: massma/propagator-bbvi

max_issues_repo_name: massma/propagator-bbvi

max_forks_repo_name: massma/propagator-bbvi


text:
module Spec_3_4_3 where

import Test.Hspec
import Numeric.LinearAlgebra

import Section_3_4_3 (initNetwork, forward)

spec :: Spec
spec = do
  describe "forward" $ do
    it "(1.0, 0.5) -> (0.326, 0.696)" $ do
        let x = vector [1.0, 0.5]
            answer = vector [0.31682708, 0.69627909]
            prediction = forward initNetwork x
            err = 0.0001
        norm_2 (answer - prediction) < err `shouldBe` True


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
module Data.Array.ShapedS.MatMul(matMul) where
import Prelude hiding ((<>))
import GHC.Stack
import GHC.TypeLits
import Data.Array.Internal(valueOf)
import Data.Array.ShapedS
import Numeric.LinearAlgebra as N

matMul :: forall m n o a .
          (HasCallStack, N.Numeric a, KnownNat m, KnownNat n, KnownNat o) =>
          Array [m, n] a -> Array [n, o] a -> Array [m, o] a
matMul x y =
  let n = valueOf @n
      o = valueOf @o
      x' = N.reshape n $ toVector x
      y' = N.reshape o $ toVector y
      xy' = x' N.<> y'
      xy = fromVector $ N.flatten xy'
  in  xy

max_stars_repo_name: augustss/orthotope-hmatrix

max_issues_repo_name: augustss/orthotope-hmatrix

max_forks_repo_name: augustss/orthotope-hmatrix


text:
{-# LANGUAGE ScopedTypeVariables, TypeFamilies, TupleSections #-}

module Fuml.Core where

import qualified Data.Vector.Storable as VS
import Numeric.LinearAlgebra
import Data.List (nub)
import Lens.Micro
import Control.Monad.Identity

data Weighted o = Weighted Double o

-- |The result of running a predictive model
data Predict p a = Predict
  { model :: p -- ^ the internal state of the model
  , predict :: Vector Double -> a -- ^ the prediction function
  }

instance Functor (Predict p) where
  fmap f (Predict m predf) = Predict m (f . predf)

newtype Supervisor m o p a
  = Supervisor { runSupervisor :: Maybe p
                               -> [(Vector Double, o)]
                               -> m (Predict p a) }

instance Functor m => Functor (Supervisor m o p) where
  fmap f (Supervisor sf) = Supervisor $ \mp d -> fmap (fmap f) $ sf mp d

-- | Helper function for running simple Supervisors in the Identity monad
runSupervisor' :: Supervisor Identity o p a -- ^ the 'Supervisor' value
               -> [(Vector Double, o)] -- ^ the dataset
               -> Predict p a -- ^ the 'Predict' value
runSupervisor' (Supervisor sf) = runIdentity . sf Nothing

oneVsRest :: (Monad m, Eq a, Functor m) => Supervisor m Bool p b -> Supervisor m a [(a,p)] [(a,b)]
oneVsRest subsuper  = Supervisor $ \_ theData -> do
  let classes = nub $ map snd theData
      boolize c (v, c1) = (v, c == c1)
      train c = fmap (c,) $ runSupervisor subsuper Nothing $ map (boolize c) theData
  models <-  mapM train classes
  return $ Predict (map (over _2 model) models) $ \v -> map (\(c,pr) -> (c, predict pr v)) models

(~~) :: (a -> o) -> [a -> Double] -> [a] -> [(Vector Double, o)]
fy ~~ fxs = let f w = (VS.fromList (map ($w) fxs) , fy w)
            in map f

--prepare :: (a -> o) -> [a -> Double] -> [a] -> [(Vector Double, o)]
--prepare fo fxs d = fo ~~ fxs $ d

max_stars_repo_name: ekalosak/open

max_issues_repo_name: ekalosak/open

max_forks_repo_name: ekalosak/open


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE NoImplicitPrelude #-}
module Main where

import qualified Prelude as P
import Data.Array.Accelerate.Data.Complex as A
import Data.Array.Accelerate as A
-- import qualified Data.Array.Accelerate.Fourier.Adhoc as A
import qualified Data.Array.Accelerate.Interpreter as I
import Data.Array.Accelerate.Interpreter as I
import Data.Array.Accelerate.Data.Fold
import Data.Array.Accelerate.Data.Monoid
import qualified Data.List as L
import qualified Data.Array.Accelerate.LLVM.PTX as PTX

import qualified Data.Csv as CSV
import qualified Data.ByteString.Lazy as BS
import qualified Statistics.Sample as S
import qualified Data.Vector as V

type Arr a = Acc (Vector a)

mean :: (Floating a, FromIntegral Int a) => Arr a -> Acc (Scalar a)
mean xs = unit $ the (fold (+) 0 xs) / (fromIntegral $ length xs)

smooth :: (Floating a, FromIntegral Int a) => Arr a -> Exp Int -> Arr a
smooth xs res = error "LOL"

sumF :: (Elt a, Lift Exp a, P.Num (Exp a)) => Fold (Exp a) (Exp a)
sumF = Fold (lift . Sum) (getSum . unlift)

lengthF :: (Elt a, Lift Exp a, P.Num (Exp a)) => Fold (Exp a) (Exp a)
lengthF = Fold (\_ -> 1) (getSum . unlift)

meanF :: (Elt a, Lift Exp a, P.Fractional (Exp a)) => Fold (Exp a) (Exp a)
meanF = sumF / lengthF

valF :: (Elt a, Monoid (Exp a)) => Fold (Exp a) (Exp a)
valF = Fold P.id P.id

-- | E(X^2)
sumOfSquaresF :: (Elt a, Lift Exp a, P.Num (Exp a)) => Fold (Exp a) (Exp a)
sumOfSquaresF = Fold (lift . Sum . (^(2 :: Exp Int))) (getSum . unlift)

-- | sqrt (E(X^2) - E(X) ^ 2)
standardDeviationF :: (Elt a, Lift Exp a, P.Floating (Exp a)) => Fold (Exp a) (Exp a)
standardDeviationF = sqrt $ (sumOfSquaresF / lengthF) - (meanF * meanF)

-- | E(X - mu) ^ 4
-- differenceQuadsF :: (Elt a, Lift Exp a, P.Fractional (Exp a), Monoid (Exp a)) => Fold (Exp a) (Exp a)
-- differenceQuadsF = (\m -> Fold (\x -> lift $ Sum ((x - m) ^ e)) (getSum . unlift)) P.<$> meanF
--   where
--     e = 4 :: Exp Int

-- | E((X - u)^2)^2
-- otherThing :: (Elt a, Lift Exp a, P.Fractional (Exp a), Monoid (Exp a)) => Fold (Exp a) (Exp a)
-- otherThing = Fold (lift . Sum . (^(2 :: Exp Int))) (getSum . unlift)

-- kurtosisF :: (Elt a, Lift Exp a, P.Fractional (Exp a), Monoid (Exp a)) => Fold (Exp a) (Exp a)
-- kurtosisF = (valF - meanF)

-- | Kurt(X) = E(X - u)^4/(E((X - u)^2)^2)
-- (mean (map (\x -> x - u) xs)) ^ four / ((mean (map (\x -> (x - u) ^ two) xs))^two)
kurtosis :: (Floating a, Elt a, FromIntegral Int a) => Arr a -> Exp a
kurtosis xs = (u4 / (variance P.^ 2))
  where
    u4 = the $ compute $ mean $ map (\x -> (x - u) ^ (4 :: Exp Int)) xs
    variance = the $ compute $ mean $ map (\x -> (x - u) ^ (2 :: Exp Int)) xs
    u = the $ compute $ mean xs

diffs :: (Elt a, P.Num (Exp a)) => Arr a -> Arr a
diffs xs = zipWith (-) (tail xs) xs

roughness :: (Elt a, P.Floating (Exp a), FromIntegral Int a, Lift Exp a) => Arr a -> Acc (Scalar a)
roughness xs = compute $ runFold standardDeviationF $ diffs xs

-- smoothSimple :: (Elt a, Floating a, FromIntegral Int a, Ord a) => Arr a -> Exp Int -> Arr a
-- smoothSimple xs resolution =
--   ifThenElse (resolution < length xs)
--     smaLargeResolution
--     (sma bestWindowSize 1 xs)
--   where
--     -- THIS IS WRONG. IT NEEDS TO MINIMIZE THE ROUGHNESS
--     bestWindowSize =  the $ minimum $ map fst $ afst $ filter acceptableWindow $ map (go xs) ws
--     acceptableWindow b =
--       (r < originalRoughness) && (k >= originalKurt)
--       where
--         r = fst $ snd b
--         k = snd $ snd b
--     smaLargeResolution = sma (length xs `div` resolution) (length xs `div` resolution) xs
--     originalKurt = kurtosis xs
--     originalRoughness = the $ roughness xs
--     end :: Exp Int
--     end = length xs `div` 10
--     ws :: Arr Int
--     ws = enumFromN (index1 $ end - 2) 2

-- smoothSimple' :: (Elt a, Floating a, FromIntegral Int a, Ord a) => Arr a -> Exp Int -> Arr a
-- smoothSimple' xs resolution = sma bestWindowSize 1 xs
--   where
--     bestWindowSize =
--       fst $ while (\a -> fst a >= 2) godude $ lift ((length xs `div` resolution), (originalRoughness, 1 :: Exp Int))
--     originalKurt = kurtosis xs
--     originalRoughness = the $ roughness xs
--     godude a =
--       let
--         w = fst a
--         minObj = fst $ snd a
--         smoothed = sma w 1 xs
--         r = the $ roughness smoothed
--         kurt = kurtosis smoothed
--       in ifThenElse (r < minObj && kurt >= originalKurt) (lift (w-1, (r, w))) (lift (w-1, snd a))

-- go :: (Elt a, P.Floating (Exp a), FromIntegral Int a) => Arr a -> Exp Int -> Exp (Int, (a, a))
-- go xs w = lift (w,(the $ roughness smoothed, kurtosis smoothed))
--   where
--     smoothed = sma w 1 xs

sma :: (Floating a, FromIntegral Int a) => Exp Int -> Exp Int -> Arr a -> Arr a
sma range slide xs = mavg range xs
  where
    skip ys =
      map snd. afst
        $ filter (\y -> fst y `mod` slide == 0)
          $ imap (\i y -> lift (unindex1 i, y)) ys

computeStuff :: (Floating a, FromIntegral Int a, Lift Exp a) => Exp Int -> Exp Int -> Arr a -> Exp (a, a)
computeStuff resolution window xs = lift (the $ roughness smoothed, kurtosis smoothed)
  where
    smoothed = sma window resolution xs

kurtAndRoughness :: (Floating a, FromIntegral Int a, Lift Exp a) => Arr a -> Exp (a, a)
kurtAndRoughness xs = lift (the $ roughness xs, kurtosis xs)

mavg :: (Floating a, FromIntegral Int a) => Exp Int -> Arr a -> Arr a
mavg k xs =
  map (/ fromIntegral k)
    $ scanl (+) (the $ sum h)
      $ zipWith (-) t xs
  where
    h = take k xs
    t = drop k xs

-- transform :: Arr (Complex Double) -> Arr (Complex Double)
-- transform = A.transform A.forward

-- inverseTransform :: Arr (Complex Double) -> Arr (Complex Double)
-- inverseTransform = A.transform A.inverse

binarySearch :: Int -> Int -> Arr (Complex Double) -> a -> a -> a -> ()
binarySearch h t xs minObj originalKurt windowSize = ()

main :: P.IO ()
main = do
  P.Right xs <- CSV.decode CSV.HasHeader P.<$> BS.readFile "data.csv" :: P.IO (P.Either P.String (V.Vector (P.String, Double)))

  let
    -- use accelerate-io
    v = V.toList $ P.snd P.<$> xs :: [Double]
    vec = use $ fromList (Z :. (V.length xs)) v
    ws = P.reverse [2..V.length xs `div` 10]
    [(minObj, originalKurt)] = toList $ PTX.runN $ unit $ kurtAndRoughness vec
    runOnWindow w = P.head $ toList $ PTX.runN $ unit $ computeStuff 1 (constant w) vec
    roughsAndKurts = (\w -> (w, runOnWindow w)) P.<$> ws
    bestWindow@(_, w) = L.foldl' goman (minObj, 1) roughsAndKurts
    goman old@(obj, _) (w,(r, k)) = if (r P.< obj P.&& k P.>= originalKurt) then (r,w) else old
    smoothed = PTX.run $ sma (constant w) 1 vec
  P.print $ L.take 100 roughsAndKurts
  P.print $ (minObj, originalKurt)
  P.print $ V.length xs
  P.print smoothed
  -- P.print $ S.kurtosis $ (V.fromList v :: V.Vector Double)
    -- bestWindowSize =  the $ minimum $ map fst $ afst $ filter acceptableWindow $ map (go xs) ws
    -- acceptableWindow b =
    --   (r < originalRoughness) && (k >= originalKurt)
    --   where
    --     r = fst $ snd b
    --     k = snd $ snd b
    -- smaLargeResolution = sma (length xs `div` resolution) (length xs `div` resolution) xs
    -- originalKurt = kurtosis xs
    -- originalRoughness = the $ roughness xs
    -- end :: Exp Int
    -- end = length xs `div` 10
    -- ws :: Arr Int
    -- ws = enumFromN (index1 $ end - 2) 2

mavg' :: P.Fractional b => Int -> [b] -> [b]
mavg' k lst = P.map (/ P.fromIntegral k) $ L.scanl' (+) (P.sum h) $ P.zipWith (-) t lst
  where (h, t) = L.splitAt k lst

max_stars_repo_name: cotrone/hs-asap

max_issues_repo_name: cotrone/hs-asap

max_forks_repo_name: cotrone/hs-asap