text:
module Render (cfileSink) where

import           Control.Monad
import           Data.ByteString.Builder
import           Data.Complex
import           IQ
import           Pipes
import qualified Pipes.Prelude as P
import           System.IO

foldChunk :: (Monad m, Monoid a) => Int -> Pipe a a m ()
foldChunk n = forever $ go n mempty
  where go 0 x = yield x
        go n x = do
          x' <- await
          go (n - 1) $ x <> x'

renderSample :: IQ -> Builder
renderSample (real :+ imag) = floatLE real <> floatLE imag

cfileSink :: MonadIO m => Handle -> Consumer IQ m ()
cfileSink h = P.map renderSample
           >-> foldChunk 100
           >-> P.mapM_ (liftIO . hPutBuilder h)
{-# INLINE cfileSink #-}

max_stars_repo_name: hexagonal-sun/ayeQ

max_issues_repo_name: hexagonal-sun/ayeQ

max_forks_repo_name: hexagonal-sun/ayeQ


text:
{-# LANGUAGE DeriveDataTypeable    #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE KindSignatures        #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PatternSynonyms       #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE RecordWildCards       #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeInType            #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE ViewPatterns          #-}

module Backprop.Learn.Model.Stochastic (
    dropout
  , rreLU
  , injectNoise, applyNoise
  , injectNoiseR, applyNoiseR
  ) where

import           Backprop.Learn.Model.Function
import           Backprop.Learn.Model.Types
import           Control.Monad.Primitive
import           Data.Bool
import           GHC.TypeNats
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop
import           Numeric.LinearAlgebra.Static.Vector
import qualified Data.Vector.Storable.Sized            as SVS
import qualified Statistics.Distribution               as Stat
import qualified System.Random.MWC                     as MWC
import qualified System.Random.MWC.Distributions       as MWC

-- | Dropout layer.  Parameterized by dropout percentage (should be between
-- 0 and 1).
--
-- 0 corresponds to no dropout, 1 corresponds to complete dropout of all
-- nodes every time.
dropout
    :: KnownNat n
    => Double
    -> Model 'Nothing 'Nothing (R n) (R n)
dropout r = Func
    { runFunc      = (auto (realToFrac (1 - r)) *)
    , runFuncStoch = \g x -> do
        (x *) . auto . vecR <$> SVS.replicateM (mask g)
    }
  where
    mask :: PrimMonad m => MWC.Gen (PrimState m) -> m Double
    mask = fmap (bool 1 0) . MWC.bernoulli r

-- | Random leaky rectified linear unit
rreLU
    :: (Stat.ContGen d, Stat.Mean d, KnownNat n)
    => d
    -> Model 'Nothing 'Nothing (R n) (R n)
rreLU d = Func
    { runFunc      = vmap' (preLU v)
    , runFuncStoch = \g x -> do
        α <- vecR <$> SVS.replicateM (Stat.genContVar d g)
        pure (zipWithVector preLU (constVar α) x)
    }
  where
    v :: BVar s Double
    v = auto (Stat.mean d)

-- | Inject random noise.  Usually used between neural network layers, or
-- at the very beginning to pre-process input.
--
-- In non-stochastic mode, this adds the mean of the distribution.
injectNoise
    :: (Stat.ContGen d, Stat.Mean d, Fractional a)
    => d
    -> Model 'Nothing 'Nothing a a
injectNoise d = Func
    { runFunc = (realToFrac (Stat.mean d) +)
    , runFuncStoch = \g x -> do
        e <- Stat.genContVar d g
        pure (realToFrac e + x)
    }

-- | 'injectNoise' lifted to 'R'
injectNoiseR
    :: (Stat.ContGen d, Stat.Mean d, KnownNat n)
    => d
    -> Model 'Nothing 'Nothing (R n) (R n)
injectNoiseR d = Func
    { runFunc = (realToFrac (Stat.mean d) +)
    , runFuncStoch  = \g x -> do
        e <- vecR <$> SVS.replicateM (Stat.genContVar d g)
        pure (constVar e + x)
    }

-- | Multply by random noise.  Can be used to implement dropout-like
-- behavior.
--
-- In non-stochastic mode, this scales by the mean of the distribution.
applyNoise
    :: (Stat.ContGen d, Stat.Mean d, Fractional a)
    => d
    -> Model 'Nothing 'Nothing a a
applyNoise d = Func
    { runFunc = (realToFrac (Stat.mean d) *)
    , runFuncStoch  = \g x -> do
        e <- Stat.genContVar d g
        pure (realToFrac e * x)
    }

-- | 'applyNoise' lifted to 'R'
applyNoiseR
    :: (Stat.ContGen d, Stat.Mean d, KnownNat n)
    => d
    -> Model 'Nothing 'Nothing (R n) (R n)
applyNoiseR d = Func
    { runFunc = (realToFrac (Stat.mean d) *)
    , runFuncStoch  = \g x -> do
        e <- vecR <$> SVS.replicateM (Stat.genContVar d g)
        pure (constVar e * x)
    }

ext: hs

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# LANGUAGE OverloadedStrings #-}

module JsonParser where

import           Data.Aeson              hiding ( Series )
import           Data.Aeson.Types               ( parseMaybe )
import qualified Data.ByteString.Lazy          as BSL
import           Data.Complex                   ( Complex((:+)) )

import           Form

data Vec = Vec {
    number :: Int,
    real   :: Double,
    imag   :: Double
} deriving (Show)

instance ToJSON Vec where
    toJSON (Vec number real imag) =
        object ["number" .= number, "real" .= real, "imag" .= imag]

instance FromJSON Vec where
    parseJSON = withObject "Vec"
        $ \v -> Vec <$> v .: "number" <*> v .: "real" <*> v .: "imag"

newtype Series = Series {
    series :: [Vec]
} deriving (Show)

instance ToJSON Series where
    toJSON (Series series) = object ["series" .= series]

instance FromJSON Series where
    parseJSON = withObject "Series" $ \s -> Series <$> s .: "series"

instance Semigroup Series where
    (Series xs) <> (Series ys) = Series (xs <> ys)

instance Monoid Series where
    mempty = Series []

convert :: Vec -> Form
convert (Vec n r i) = form (r :+ i) n

getForm :: Series -> [Form]
getForm (Series xs) = map convert xs

jsonToCoord :: FilePath -> IO [(Int, Int)]
jsonToCoord path = formsToCoords <$> jsonToForm path

jsonToForm :: FilePath -> IO [Form]
jsonToForm path = do
    json <- BSL.readFile path
    case decode json of
        Just js -> return $ getForm js
        Nothing -> error $ "invalid json file: " ++ path

jsonToCoords :: [FilePath] -> IO [(Int, Int)]
jsonToCoords = foldr (mappend . jsonToCoord) mempty


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE Strict           #-}
module FourierMethod.FourierSeries2D
  ( module FourierMethod.FourierSeries2D
  , module FourierMethod.BlockCudaMatrix
  , module FourierMethod.BlockMatrixAcc
  ) where

import           Control.DeepSeq
import           Control.Monad                               as M
import           Control.Monad.IO.Class
import           Control.Monad.Trans.Resource
import qualified Data.Array.Accelerate                       as A
import           Data.Array.Accelerate.LLVM.PTX              as A
import           Data.Array.Accelerate.Numeric.LinearAlgebra as A
import           Data.Array.Repa                             as R
import           Data.Complex
import           Data.Conduit                                as C
import           Data.Conduit.List                           as CL
import           Data.List                                   as L
import           Data.Vector.Storable                        as VS
import           Data.Vector.Unboxed                         as VU
import           Foreign.CUDA.Driver                         as CUDA
import           FourierMethod.BlockCudaMatrix
import           FourierMethod.BlockMatrixAcc
import           FourierMethod.FourierSeries2DAcc
import           Utils.BLAS
import           Utils.Distribution
import           Utils.List
import           Utils.Parallel                              hiding ((.|))
import           Utils.Time

{-# INLINE harmonicMatPTX #-}
harmonicMatPTX ::
     ( Storable e
     , A.Elt (Complex e)
     , A.Elt e
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     )
  => Int
  -> e
  -> e
  -> e
  -> PTX
  -> [(Int, Int)]
  -> CuMat (Complex e)
harmonicMatPTX numPoints period delta deltaFreq ptx r2Freqs =
  let rows = L.length r2Freqs
  in CuMat rows (numPoints ^ 2) .
     CuVecHost .
     VS.fromList .
     A.toList .
     run1With ptx (harmonicAcc numPoints period delta deltaFreq) .
     A.fromList (A.Z A.:. (L.length r2Freqs)) $
     r2Freqs

-- create column-major inverse harmonics
{-# INLINE inverseHarmonicMatPTX #-}
inverseHarmonicMatPTX ::
     ( Storable e
     , A.Elt (Complex e)
     , A.Elt e
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     )
  => Int
  -> e
  -> e
  -> PTX
  -> [(Int, Int)]
  -> CuMat (Complex e)
inverseHarmonicMatPTX numFreqs period delta ptx r2Positions =
  let cols = L.length r2Positions
  in CuMat (numFreqs ^ 2) cols .
     CuVecHost .
     VS.fromList .
     A.toList .
     run1With ptx (inverseHarmonicAcc numFreqs period delta) .
     A.fromList (A.Z A.:. (L.length r2Positions)) $
     r2Positions

{-# INLINE inverseHarmonicMatPTX1 #-}
inverseHarmonicMatPTX1 ::
     ( Storable e
     , A.Elt (Complex e)
     , A.Elt e
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     )
  => Int
  -> e
  -> e
  -> PTX
  -> [(Int, Int)]
  -> CuMat (Complex e)
inverseHarmonicMatPTX1 numFreqs period delta ptx r2Positions =
  let rows = L.length r2Positions
  in CuMat rows (numFreqs ^ 2) .
     CuVecHost .
     VS.fromList .
     A.toList .
     run1With ptx (inverseHarmonicAcc1 numFreqs period delta) .
     A.fromList (A.Z A.:. (L.length r2Positions)) $
     r2Positions

createHarmonicMatriesGPU ::
     ( Storable e
     , A.Elt (Complex e)
     , A.Elt e
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     )
  => [PTX]
  -> Int
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> IO [[CuMat (Complex e)]]
createHarmonicMatriesGPU ptxs numBatch numPoints numFreqs period delta deltaFreq = do
  let idxs =
        L.map (divideListN numBatch) . divideListN (L.length ptxs) $
        [ (xFreq, yFreq)
        | xFreq <- getListFromNumber numFreqs
        , yFreq <- getListFromNumber numFreqs
        ]
  let output =
        parZipWith
          rdeepseq
          (\ptx -> L.map (harmonicMatPTX numPoints period delta deltaFreq ptx))
          ptxs $
        idxs
  return output

createInverseHarmonicMatriesGPU ::
     ( Storable e
     , A.Elt (Complex e)
     , A.Elt e
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     )
  => [PTX]
  -> Int
  -> Int
  -> Int
  -> e
  -> e
  -> IO [[CuMat (Complex e)]]
createInverseHarmonicMatriesGPU ptxs numBatch numPoints numFreqs period delta = do
  let idxs =
        L.map (divideListN numBatch) . divideListN (L.length ptxs) $
        [ (x, y)
        | x <- getListFromNumber numPoints
        , y <- getListFromNumber numPoints
        ]
  return .
    parZipWith
      rdeepseq
      (\ptx -> L.map (inverseHarmonicMatPTX numFreqs period delta ptx))
      ptxs $
    idxs

computeFourierCoefficientsR2 ::
     ( Storable e
     , CUBLAS (Complex e)
     , Unbox e
     , RealFloat e
     , A.Elt e
     , A.Elt (Complex e)
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     )
  => [Int]
  -> [PTX]
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> Int
  -> [CuMat (Complex e)]
  -> IO (R.Array U DIM3 (Complex e))
computeFourierCoefficientsR2 deviceIDs ptxs numFreqs numPoints period delta deltaFreq numBatchR2Freqs xs = do
  let simpsonNorm = (delta / 3) ^ 2 :+ 0
      cols = L.foldl' (\s mat -> s + getColsCuMat mat) 0 xs
  harmonics <-
    createHarmonicMatriesGPU
      ptxs
      numBatchR2Freqs
      numPoints
      numFreqs
      period
      delta
      deltaFreq
  coefs <- blockMatrixMultiply2 True deviceIDs harmonics xs
  return .
    fromUnboxed (Z :. cols :. numFreqs :. numFreqs) .
    VU.map (* simpsonNorm) . VS.convert . getHostCuMat $
    coefs

computeFourierSeriesR2 ::
     ( Storable e
     , CUBLAS (Complex e)
     , Unbox e
     , RealFloat e
     , A.Elt e
     , A.Elt (Complex e)
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     )
  => [Int]
  -> Int
  -> Int
  -> e
  -> [[CuMat (Complex e)]]
  -> [CuMat (Complex e)]
  -> IO (R.Array U DIM3 (Complex e))
computeFourierSeriesR2 deviceIDs numFreqs numPoints period inverseHarmonics xs = do
  let rows = L.foldl' (\s mat -> s + getRowsCuMat mat) 0 xs
  series <- blockMatrixMultiply1 False deviceIDs xs inverseHarmonics
  return .
    fromUnboxed (Z :. rows :. numPoints :. numPoints) .
    VS.convert . getHostCuMat $
    series

{-# INLINE applyGaussian #-}
applyGaussian ::
     (R.Source s (Complex e), Unbox e, RealFloat e)
  => e
  -> e
  -> R.Array s DIM3 (Complex e)
  -> R.Array D DIM3 (Complex e)
applyGaussian deltaFreq std arr =
  let (Z :. _ :. numFreq :. _) = extent arr
      freq = div numFreq 2
      gaussianCoefficients =
        computeUnboxedS . R.fromFunction (Z :. numFreq :. numFreq) $ \(Z :. xFreq :. yFreq) ->
          gaussian2D
            (deltaFreq * (fromIntegral $ xFreq - freq))
            (deltaFreq * (fromIntegral $ yFreq - freq))
            std :+
          0
  in R.traverse2 arr gaussianCoefficients const $ \fC fG idx@(Z :. _ :. i :. j) ->
       fC idx * fG (Z :. i :. j)

-- Stream
{-# INLINE indexSource #-}
indexSource :: Int -> Int -> ConduitT () [(Int, Int)] (ResourceT IO) ()
indexSource numIndex numBatch =
  let xs =
        divideListN numBatch $
        [ (xFreq, yFreq)
        | xFreq <- getListFromNumber numIndex
        , yFreq <- getListFromNumber numIndex
        ]
  in CL.sourceList xs

{-# INLINE fourierCoefficientsConduit #-}
fourierCoefficientsConduit ::
     ( Storable e
     , CUBLAS (Complex e)
     , Unbox e
     , RealFloat e
     , A.Elt e
     , A.Elt (Complex e)
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     , Show e
     )
  => [Int]
  -> [PTX]
  -> Int
  -> e
  -> e
  -> e
  -> [CuMat (Complex e)]
  -> ConduitT [(Int, Int)] (CuMat (Complex e)) (ResourceT IO) ()
fourierCoefficientsConduit deviceIDs ptxs numPoints period delta deltaFreq xs =
  awaitForever $ \idx' -> do
    let idxs = divideListN (L.length deviceIDs) idx'
        harmonics =
          parMap
            rdeepseq
            (uncurry (harmonicMatPTX numPoints period delta deltaFreq)) $
          L.zip ptxs idxs
    liftIO $ printCurrentTime ""
    coefs <- liftIO $ blockMatrixMultiply3 False deviceIDs harmonics xs
    yield $!! coefs

{-# INLINE fourierCoefficientsSink #-}
fourierCoefficientsSink ::
     (Storable e, Unbox e, RealFloat e, Show e)
  => Int
  -> e
  -> ConduitT (CuMat (Complex e)) Void (ResourceT IO) ((R.Array U DIM3 (Complex e)))
fourierCoefficientsSink numR2Freqs delta = do
  xs <- CL.consume
  let cols = getColsCuMat . L.head $ xs
      simpsonNorm = (delta / 3) ^ 2 :+ 0
  return .
    fromUnboxed (Z :. cols :. numR2Freqs :. numR2Freqs) .
    VU.map (* simpsonNorm) .
    VS.convert . getHostCuMat . transposeCuMat . concatCuMat $
    xs

computeFourierCoefficientsR2Stream ::
     ( Storable e
     , CUBLAS (Complex e)
     , Unbox e
     , RealFloat e
     , A.Elt e
     , A.Elt (Complex e)
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     )
  => [Int]
  -> [PTX]
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> Int
  -> [CuMat (Complex e)]
  -> IO (R.Array U DIM3 (Complex e))
computeFourierCoefficientsR2Stream deviceIDs ptxs numFreqs numPoints period delta deltaFreq numBatch xs =
  runConduitRes $
  indexSource numFreqs numBatch .|
  fourierCoefficientsConduit deviceIDs ptxs numPoints period delta deltaFreq xs .|
  fourierCoefficientsSink numFreqs delta

{-# INLINE fourierSeriesConduit #-}
fourierSeriesConduit ::
     ( Storable e
     , CUBLAS (Complex e)
     , Unbox e
     , RealFloat e
     , A.Elt e
     , A.Elt (Complex e)
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     , Show e
     )
  => [Int]
  -> [PTX]
  -> Int
  -> e
  -> e
  -> [CuMat (Complex e)]
  -> ConduitT [(Int, Int)] (CuMat (Complex e)) (ResourceT IO) ()
fourierSeriesConduit deviceIDs ptxs numFreqs period delta xs =
  awaitForever $ \idx' -> do
    let idxs = divideListN (L.length deviceIDs) idx'
        harmonics =
          parMap
            rdeepseq
            (uncurry (inverseHarmonicMatPTX1 numFreqs period delta)) $
          L.zip ptxs idxs
    liftIO $ printCurrentTime "fourierSeriesConduit"
    coefs <- liftIO $ blockMatrixMultiply3 False deviceIDs harmonics xs
    yield $!! coefs

{-# INLINE fourierSeriesSink #-}
fourierSeriesSink ::
     (Storable e, Unbox e, RealFloat e, Show e)
  => Int
  -> ConduitT (CuMat (Complex e)) Void (ResourceT IO) ((R.Array U DIM3 (Complex e)))
fourierSeriesSink numPoints = do
  xs <- CL.consume
  let cols = getColsCuMat . L.head $ xs
  return .
    fromUnboxed (Z :. cols :. numPoints :. numPoints) .
    VS.convert . getHostCuMat . transposeCuMat . concatCuMat $
    xs

computeFourierSeriesR2Stream ::
     ( Storable e
     , CUBLAS (Complex e)
     , Unbox e
     , RealFloat e
     , A.Elt e
     , A.Elt (Complex e)
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     )
  => [Int]
  -> [PTX]
  -> Int
  -> Int
  -> e
  -> e
  -> Int
  -> [CuMat (Complex e)]
  -> IO (R.Array U DIM3 (Complex e))
computeFourierSeriesR2Stream deviceIDs ptxs numFreqs numPoints period delta numBatch xs =
  runConduitRes $
    indexSource numPoints numBatch .|
    fourierSeriesConduit deviceIDs ptxs numFreqs period delta xs .|
    fourierSeriesSink numPoints

{-# INLINE fourierSeriesConduitAcc #-}
fourierSeriesConduitAcc ::
     ( A.Floating e
     , A.Elt (Complex e)
     , A.FromIntegral Int e
     , Numeric (Complex e)
     , Unbox e
     )
  => [PTX]
  -> Int
  -> e
  -> e
  -> A.Acc (Matrix (Complex e))
  -> ConduitT [(Int, Int)] (VU.Vector (Complex e)) (ResourceT IO) ()
fourierSeriesConduitAcc ptxs numFreqs period delta x =
  awaitForever $ \idx' -> do
    let idxs = divideListN (L.length ptxs) idx'
        harmonics =
          L.map
            (\r2Positions ->
               inverseHarmonicAcc2 numFreqs period delta .
               A.fromList (A.Z A.:. L.length r2Positions) $
               r2Positions) $
          idxs
    liftIO $ printCurrentTime "fourierSeriesConduit"
    yield $!! blockMatrixMultiply ptxs harmonics x

{-# INLINE fourierSeriesSinkAcc #-}
fourierSeriesSinkAcc ::
     (Unbox e)
  => Int
  -> Int
  -> ConduitT (VU.Vector (Complex e)) Void (ResourceT IO) (R.Array U DIM3 (Complex e))
fourierSeriesSinkAcc numPoints cols = do
  xs <- CL.consume
  computeP .
    R.backpermute
      (Z :. cols :. numPoints :. numPoints)
      (\(Z :. a :. b :. c) -> Z :. b :. c :. a) .
    fromUnboxed (Z :. numPoints :. numPoints :. cols) . VU.concat $
    xs

computeFourierSeriesR2StreamAcc ::
     ( A.Floating e
     , A.Elt (Complex e)
     , A.FromIntegral Int e
     , Numeric (Complex e)
     , Unbox e
     )
  => [PTX]
  -> Int
  -> Int
  -> Int
  -> e
  -> e
  -> Int
  -> Acc (Matrix (Complex e))
  -> IO (R.Array U DIM3 (Complex e))
computeFourierSeriesR2StreamAcc ptxs numFreqs numPoints cols period delta numBatch x =
  runConduitRes $
  indexSource numPoints numBatch .|
  fourierSeriesConduitAcc ptxs numFreqs period delta (A.compute x) .|
  fourierSeriesSinkAcc numPoints cols

{-# INLINE fourierSeriesConduitAcc' #-}
fourierSeriesConduitAcc' ::
     ( A.Floating e
     , A.Elt (Complex e)
     , A.FromIntegral Int e
     , Numeric (Complex e)
     , Unbox e
     , Prelude.Num e
     )
  => [PTX]
  -> Int
  -> e
  -> e
  -> A.Acc (Matrix (Complex e))
  -> ConduitT [(Int, Int)] (VU.Vector (Complex e)) (ResourceT IO) ()
fourierSeriesConduitAcc' ptxs numFreqs period delta x =
  awaitForever $ \idx' -> do
    let idxs = divideListN (L.length ptxs) idx'
        harmonics =
          L.map
            (\r2Positions ->
               inverseHarmonicAcc2' numFreqs period delta .
               A.fromList (A.Z A.:. L.length r2Positions) $
               r2Positions) $
          idxs
    liftIO $ printCurrentTime "fourierSeriesConduit"
    yield $ blockMatrixMultiply ptxs harmonics x

computeFourierSeriesR2StreamAcc' ::
     ( A.Floating e
     , A.Elt (Complex e)
     , A.FromIntegral Int e
     , Numeric (Complex e)
     , Unbox e
     , Prelude.Num e
     )
  => [PTX]
  -> Int
  -> Int
  -> Int
  -> e
  -> e
  -> Int
  -> Acc (Matrix (Complex e))
  -> IO (R.Array U DIM3 (Complex e))
computeFourierSeriesR2StreamAcc' ptxs numFreqs numPoints cols period delta numBatch x =
  runConduitRes $
  indexSource numPoints numBatch .|
  fourierSeriesConduitAcc' ptxs numFreqs period delta (A.compute x) .|
  fourierSeriesSinkAcc numPoints cols

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:
module Math.FROG.Retrieval
( retrieve
, calcLoss
) where

import qualified Numeric.LinearAlgebra as LA
import qualified Math.Optimization.SPSA as SPSA

import qualified Math.FROG.Tools as FR
import Math.FROG.Types

retrieve :: Nonlinearity -> Trace -> IO (LA.Vector Double)
retrieve nonlin measurement = do
    let niter = 50000
    let a = 0.2
    let c = 0.00001

    let sz = LA.rows measurement
    let gainA = SPSA.standardAk a (niter `quot` 10) 0.602
    let gainC = SPSA.standardCk c 0.101
    let lossFn = calcLoss nonlin measurement

    spsa <- SPSA.mkUnconstrainedSPSA lossFn gainA gainC (2*sz)
    initialGuess <- FR.mkInitialGuess sz

    return $ SPSA.optimize spsa niter initialGuess

calcLoss :: Nonlinearity -> Trace -> LA.Vector Double -> Double
calcLoss nonlin meas v = (/ (fromIntegral $ LA.rows meas)) . sqrt
                                                           . LA.sumElements $ sqDiff
    where
    field = FR.flatPolar2Complex v
    (signal, gate) = FR.mkSigGate nonlin field field
    sqDiff = (^2) . abs $ (meas - FR.mkTrace signal gate)

max_stars_repo_name: leroix/haskell-spsa-frog

max_issues_repo_name: leroix/haskell-spsa-frog

max_forks_repo_name: leroix/haskell-spsa-frog


text:
{-# LANGUAGE CApiFFI #-}

-- {-# LANGUAGE QuasiQuotes #-}
-- {-# LANGUAGE TemplateHaskell #-}

module LatticeSymmetries
  ( ls_hs_hdf5_create_dataset_u64,
    ls_hs_hdf5_create_dataset_f32,
    ls_hs_hdf5_create_dataset_f64,
    ls_hs_hdf5_create_dataset_c64,
    ls_hs_hdf5_create_dataset_c128,
    ls_hs_hdf5_write_chunk_u64,
    ls_hs_hdf5_write_chunk_f32,
    ls_hs_hdf5_write_chunk_f64,
    ls_hs_hdf5_write_chunk_c64,
    ls_hs_hdf5_write_chunk_c128,
    ls_hs_hdf5_read_chunk_u64,
    ls_hs_hdf5_read_chunk_f32,
    ls_hs_hdf5_read_chunk_f64,
    ls_hs_hdf5_read_chunk_c64,
    ls_hs_hdf5_read_chunk_c128,
    ls_hs_hdf5_get_dataset_rank,
    ls_hs_hdf5_get_dataset_shape,
  )
  where

import Data.Complex
import qualified Data.HDF5 as H5
import qualified Data.HDF5.Types as H5
import qualified Data.Vector.Storable as V
import Foreign.C.String (CString, peekCString)
import Foreign.C.Types (CUInt (..))
import Foreign.ForeignPtr
import Foreign.Marshal.Array (peekArray, pokeArray)
import Foreign.Ptr (Ptr)
import Foreign.Storable (Storable (..))

  {-
-- | Retrieve textual representation of an error
getErrorMessage ::
  -- | Error code returned by @lattice_symmetries@ library
  Int ->
  -- | Explanation of the error
  IO Text
getErrorMessage c = bracket (ls_error_to_string (fromIntegral c)) ls_destroy_string $ \s ->
  toText <$> peekCString s
-}

  {-
-- | Check the status code returned by @lattice_symmetries@ library. If it
-- indicates an error, 'LatticeSymmetriesException' is thrown.
checkStatus :: (MonadIO m, Integral a) => a -> m ()
checkStatus c
  | c == 0 = return ()
  | otherwise = do
    print =<< liftIO (getErrorMessage c')
    liftIO $ throwIO . LatticeSymmetriesException c' =<< getErrorMessage c'
  where
    c' = fromIntegral c
-}

-- loadRawBasis :: MonadIO m => Text -> m (Ptr Cspin_basis)
-- loadRawBasis path = do
--   r <- liftIO $ decodeFileWithWarnings (toString path)
--   case r of
--     Left e -> throwIO e
--     Right (warnings, (WrappedBasisSpec basisSpec)) -> do
--       mapM_ print warnings
--       toRawBasis basisSpec

-- ls_hs_internal_load_basis_from_yaml :: CString -> IO (Ptr Cspin_basis)
-- ls_hs_internal_load_basis_from_yaml c_str = loadRawBasis . toText =<< peekCString c_str
--
-- foreign export ccall "ls_hs_internal_load_basis_from_yaml"
--   ls_hs_internal_load_basis_from_yaml :: CString -> IO (Ptr Cspin_basis)

{-
isPermutation :: Integral a => NonEmpty a -> Bool
isPermutation xs = and . toList $ NonEmpty.zipWith (==) (NonEmpty.sort xs) (fromList [0 ..])

permuteArray :: Integral a => NonEmpty a -> U.Vector Int -> U.Vector Int
permuteArray permutation xs = (xs U.!!) undefined

computePeriodicity :: Integral a => NonEmpty a -> Int
computePeriodicity = undefined

checkSymmetrySpec :: SymmetrySpec -> Either Text SymmetrySpec
checkSymmetrySpec (SymmetrySpec permutation sector) = do
  when (not . isPermutation $ permutation) $
    Left $ show permutation <> " is not a permutation of " <> show [0 .. length permutation]
  undefined
-}

createDataset ::
  forall a.
  (H5.KnownDatatype a) =>
  CString ->
  CString ->
  CUInt ->
  Ptr Word64 ->
  IO ()
createDataset _fileName _datasetName _dim shapePtr = do
  fileName <- fromString <$> peekCString _fileName
  datasetName <- fromString <$> peekCString _datasetName
  shape <- fmap fromIntegral <$> peekArray (fromIntegral _dim) shapePtr
  H5.withFile fileName H5.WriteAppend $ \file -> do
    exists <- H5.exists file datasetName
    when exists $ H5.delete file datasetName
    void . join $ H5.createEmptyDataset file datasetName <$> H5.ofType @a <*> H5.ofShape shape

ls_hs_hdf5_create_dataset_u64 :: CString -> CString -> CUInt -> Ptr Word64 -> IO ()
ls_hs_hdf5_create_dataset_u64 = createDataset @Word64

ls_hs_hdf5_create_dataset_f32 :: CString -> CString -> CUInt -> Ptr Word64 -> IO ()
ls_hs_hdf5_create_dataset_f32 = createDataset @Double

ls_hs_hdf5_create_dataset_f64 :: CString -> CString -> CUInt -> Ptr Word64 -> IO ()
ls_hs_hdf5_create_dataset_f64 = createDataset @Double

ls_hs_hdf5_create_dataset_c64 :: CString -> CString -> CUInt -> Ptr Word64 -> IO ()
ls_hs_hdf5_create_dataset_c64 = createDataset @(Complex Float)

ls_hs_hdf5_create_dataset_c128 :: CString -> CString -> CUInt -> Ptr Word64 -> IO ()
ls_hs_hdf5_create_dataset_c128 = createDataset @(Complex Double)

-- foreign export ccall "ls_hs_hdf5_create_dataset_u64"
--   ls_hs_hdf5_create_dataset_u64 :: CString -> CString -> CUInt -> Ptr Word64 -> IO ()

-- foreign export ccall "ls_hs_hdf5_create_dataset_f32"
--   ls_hs_hdf5_create_dataset_f32 :: CString -> CString -> CUInt -> Ptr Word64 -> IO ()

-- foreign export ccall "ls_hs_hdf5_create_dataset_f64"
--   ls_hs_hdf5_create_dataset_f64 :: CString -> CString -> CUInt -> Ptr Word64 -> IO ()

-- foreign export ccall "ls_hs_hdf5_create_dataset_c64"
--   ls_hs_hdf5_create_dataset_c64 :: CString -> CString -> CUInt -> Ptr Word64 -> IO ()

-- foreign export ccall "ls_hs_hdf5_create_dataset_c128"
--   ls_hs_hdf5_create_dataset_c128 :: CString -> CString -> CUInt -> Ptr Word64 -> IO ()

writeDatasetChunk ::
  forall a.
  (Storable a, H5.KnownDatatype a) =>
  CString ->
  CString ->
  CUInt ->
  Ptr Word64 ->
  Ptr Word64 ->
  Ptr a ->
  IO ()
writeDatasetChunk _fileName _datasetName _dim _offsetPtr _shapePtr dataPtr = do
  fileName <- fromString <$> peekCString _fileName
  datasetName <- fromString <$> peekCString _datasetName
  let rank = fromIntegral _dim
  offset <- fromList . fmap fromIntegral <$> peekArray rank _offsetPtr
  shape <- fromList . fmap fromIntegral <$> peekArray rank _shapePtr
  let hyperslab = H5.Hyperslab offset (V.replicate rank 1) shape (V.replicate rank 1)
  fp <- newForeignPtr_ dataPtr
  let size = V.product shape
      buffer = V.unsafeFromForeignPtr0 fp size
  H5.withFile fileName H5.WriteAppend $ \file ->
    H5.open file datasetName
      >>= return
        . H5.sliceWithHyperslab hyperslab
      >>= H5.writeSelected buffer

ls_hs_hdf5_write_chunk_u64 ::
  CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr Word64 -> IO ()
ls_hs_hdf5_write_chunk_u64 = writeDatasetChunk @Word64

ls_hs_hdf5_write_chunk_f32 ::
  CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr Float -> IO ()
ls_hs_hdf5_write_chunk_f32 = writeDatasetChunk @Float

ls_hs_hdf5_write_chunk_f64 ::
  CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr Double -> IO ()
ls_hs_hdf5_write_chunk_f64 = writeDatasetChunk @Double

ls_hs_hdf5_write_chunk_c64 ::
  CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr (Complex Float) -> IO ()
ls_hs_hdf5_write_chunk_c64 = writeDatasetChunk @(Complex Float)

ls_hs_hdf5_write_chunk_c128 ::
  CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr (Complex Double) -> IO ()
ls_hs_hdf5_write_chunk_c128 = writeDatasetChunk @(Complex Double)

-- foreign export ccall "ls_hs_hdf5_write_chunk_u64"
--   ls_hs_hdf5_write_chunk_u64 :: CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr Word64 -> IO ()

-- foreign export ccall "ls_hs_hdf5_write_chunk_f64"
--   ls_hs_hdf5_write_chunk_f64 :: CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr Double -> IO ()

readDatasetChunk ::
  forall a.
  (Storable a, H5.KnownDatatype a) =>
  CString ->
  CString ->
  CUInt ->
  Ptr Word64 ->
  Ptr Word64 ->
  Ptr a ->
  IO ()
readDatasetChunk _fileName _datasetName _dim _offsetPtr _shapePtr dataPtr = do
  fileName <- fromString <$> peekCString _fileName
  datasetName <- fromString <$> peekCString _datasetName
  let rank = fromIntegral _dim
  offset <- fromList . fmap fromIntegral <$> peekArray rank _offsetPtr
  shape <- fromList . fmap fromIntegral <$> peekArray rank _shapePtr
  let hyperslab = H5.Hyperslab offset (V.replicate rank 1) shape (V.replicate rank 1)
  fp <- newForeignPtr_ dataPtr
  let size = V.product shape
      buffer = V.unsafeFromForeignPtr0 fp size
  H5.withFile fileName H5.WriteAppend $ \file ->
    H5.open file datasetName
      >>= return
        . H5.sliceWithHyperslab hyperslab
      >>= H5.readSelectedInplace buffer

ls_hs_hdf5_read_chunk_u64 ::
  CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr Word64 -> IO ()
ls_hs_hdf5_read_chunk_u64 = readDatasetChunk @Word64

ls_hs_hdf5_read_chunk_f32 ::
  CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr Float -> IO ()
ls_hs_hdf5_read_chunk_f32 = readDatasetChunk @Float

ls_hs_hdf5_read_chunk_f64 ::
  CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr Double -> IO ()
ls_hs_hdf5_read_chunk_f64 = readDatasetChunk @Double

ls_hs_hdf5_read_chunk_c64 ::
  CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr (Complex Float) -> IO ()
ls_hs_hdf5_read_chunk_c64 = readDatasetChunk @(Complex Float)

ls_hs_hdf5_read_chunk_c128 ::
  CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr (Complex Double) -> IO ()
ls_hs_hdf5_read_chunk_c128 = readDatasetChunk @(Complex Double)

-- foreign export ccall "ls_hs_hdf5_read_chunk_u64"
--   ls_hs_hdf5_read_chunk_u64 :: CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr Word64 -> IO ()

-- foreign export ccall "ls_hs_hdf5_read_chunk_f64"
--   ls_hs_hdf5_read_chunk_f64 :: CString -> CString -> CUInt -> Ptr Word64 -> Ptr Word64 -> Ptr Double -> IO ()

ls_hs_hdf5_get_dataset_rank :: CString -> CString -> IO CUInt
ls_hs_hdf5_get_dataset_rank _fileName _datasetName = do
  fileName <- fromString <$> peekCString _fileName
  datasetName <- fromString <$> peekCString _datasetName
  rank <- H5.withFile fileName H5.ReadOnly $ \file ->
    H5.open file datasetName >>= return . H5.datasetRank
  pure (fromIntegral rank)

-- foreign export ccall "ls_hs_hdf5_get_dataset_rank"
--   ls_hs_hdf5_get_dataset_rank :: CString -> CString -> IO CUInt

ls_hs_hdf5_get_dataset_shape :: CString -> CString -> Ptr Word64 -> IO ()
ls_hs_hdf5_get_dataset_shape _fileName _datasetName sizePtr = do
  fileName <- fromString <$> peekCString _fileName
  datasetName <- fromString <$> peekCString _datasetName
  shape <- H5.withFile fileName H5.ReadOnly $ \file ->
    H5.open file datasetName >>= return . H5.datasetShape
  pokeArray sizePtr (fromIntegral <$> shape)

-- foreign export ccall "ls_hs_hdf5_get_dataset_shape"
--   ls_hs_hdf5_get_dataset_shape :: CString -> CString -> Ptr Word64 -> IO ()

-- writeDatasetChunk _fileName _datasetName offset size ptr = do
--   fileName <- fromString <$> peekCString _fileName
--   datasetName <- fromString <$> peekCString _datasetName
--   H5.withFile filename H5.AppendMode $ \file ->
--     H5.open

{-
loadBasisAndHamiltonianFromYAML :: Text -> IO (SpinBasis, Operator)
loadBasisAndHamiltonianFromYAML path = do
  r <- decodeFileWithWarnings (toString path)
  case r of
    Left e -> throwIO e
    Right (warnings, (ConfigSpec basisSpec operatorSpec)) -> do
      mapM_ print warnings
      let basis = mkBasis basisSpec
      return (basis, mkOperator basis operatorSpec)
-}

{-
ls_hs_basis_and_hamiltonian_from_yaml :: CString -> Ptr SpinBasisWrapper -> Ptr OperatorWrapper -> IO ()
ls_hs_basis_and_hamiltonian_from_yaml path basisPtr operatorPtr = do
  (basis, hamiltonian) <- loadBasisAndHamiltonianFromYAML . toText =<< peekCString path
  case basis of
    SpinBasis fp -> withForeignPtr fp $ \rawPtr ->
      poke basisPtr =<< SpinBasisWrapper rawPtr <$> newStablePtr basis
  case hamiltonian of
    Operator fp -> withForeignPtr fp $ \rawPtr ->
      poke operatorPtr =<< OperatorWrapper rawPtr <$> newStablePtr hamiltonian
-}

-- void ls_hs_basis_and_hamiltonian_from_yaml(char const *path,
--                                            ls_hs_spin_basis_v1 *basis,
--                                            ls_hs_operator_v1 *hamiltonian);
-- foreign export ccall "ls_hs_basis_and_hamiltonian_from_yaml"
--   ls_hs_basis_and_hamiltonian_from_yaml :: CString -> Ptr SpinBasisWrapper -> Ptr OperatorWrapper -> IO ()

{-
ls_hs_destroy_spin_basis :: Ptr SpinBasisWrapper -> IO ()
ls_hs_destroy_spin_basis ptr = do
  peek ptr >>= \(SpinBasisWrapper _ stablePtr) -> freeStablePtr stablePtr
  poke ptr $ SpinBasisWrapper nullPtr (castPtrToStablePtr nullPtr)
-}

-- void ls_hs_destroy_spin_basis(ls_hs_spin_basis_v1 *basis);
-- foreign export ccall "ls_hs_destroy_spin_basis"
--   ls_hs_destroy_spin_basis :: Ptr SpinBasisWrapper -> IO ()

{-
ls_hs_destroy_operator :: Ptr OperatorWrapper -> IO ()
ls_hs_destroy_operator ptr = do
  peek ptr >>= \(OperatorWrapper _ stablePtr) -> freeStablePtr stablePtr
  poke ptr $ OperatorWrapper nullPtr (castPtrToStablePtr nullPtr)
-}

-- void ls_hs_destroy_operator(ls_hs_operator_v1 *op);
-- foreign export ccall "ls_hs_destroy_operator"
--   ls_hs_destroy_operator :: Ptr OperatorWrapper -> IO ()

{-
toSymmetry :: SymmetrySpec -> Symmetry
toSymmetry (SymmetrySpec p s) = mkSymmetry p s

mkSymmetry ::
  -- | Permutation
  NonEmpty Int ->
  -- | Symmetry sector
  Int ->
  Symmetry
mkSymmetry !permutation !sector
  | any (< 0) permutation =
    error $
      "invalid permutation: " <> show permutation <> "; indices must be non-negative"
  | sector < 0 =
    error $
      "invalid sector: " <> show sector <> "; expected a non-negative number"
  | otherwise = System.IO.Unsafe.unsafePerformIO $ do
    ptr <- mkObject $ \ptrPtr ->
      withArrayLen (fromIntegral <$> toList permutation) $ \n permutationPtr ->
        ls_create_symmetry ptrPtr (fromIntegral n) permutationPtr (fromIntegral sector)
    Symmetry <$> newForeignPtr ls_destroy_symmetry ptr
-}

{-
mkSymmetryGroup ::
  -- | Symmetry generators
  [Symmetry] ->
  SymmetryGroup
mkSymmetryGroup symmetries = System.IO.Unsafe.unsafePerformIO $ do
  ptr <- liftIO . mkObject $ \ptrPtr ->
    withSymmetries symmetries $ \n symmetriesPtr ->
      ls_create_group ptrPtr (fromIntegral n) symmetriesPtr
  fmap SymmetryGroup . liftIO $ newForeignPtr ls_destroy_group ptr

mkRawBasis ::
  (MonadIO m, MonadThrow m) =>
  -- | Symmetry group
  SymmetryGroup ->
  -- | Number of spins
  Int ->
  -- | Hamming weight
  Maybe Int ->
  -- | Spin inversion
  Maybe Int ->
  m SpinBasis
mkRawBasis !(SymmetryGroup g) !numberSpins !hammingWeight !spinInversion = do
  when (numberSpins <= 0) . throwIO . SpinEDException $
    "invalid number of spins: " <> show numberSpins <> "; expected a positive number"
  hammingWeight' <- case hammingWeight of
    Just x -> do
      when (x < 0) . throwIO . SpinEDException $
        "invalid Hamming weight: " <> show x <> "; expected a non-negative number"
      return $ fromIntegral x
    Nothing -> return (-1)
  spinInversion' <- case spinInversion of
    Just x -> do
      when (x /= 1 && x /= -1) . throwIO . SpinEDException $
        "invalid value for spin inversion: " <> show x <> "; expected either -1 or +1"
      return $ fromIntegral x
    Nothing -> return 0
  ptr <- liftIO . mkObject $ \ptrPtr ->
    withForeignPtr g $ \groupPtr ->
      ls_create_spin_basis ptrPtr groupPtr (fromIntegral numberSpins) hammingWeight' spinInversion'
  SpinBasis <$> liftIO (newForeignPtr ls_destroy_spin_basis ptr)
-}

{-
mkFlatBasis :: BasisSpec -> FlatSpinBasis
mkFlatBasis spec = System.IO.Unsafe.unsafePerformIO $ do
  flatBasisPtr <- alloca $ \basisPtrPtr ->
    withForeignPtr normalBasis $ \normalBasisPtr -> do
      checkStatus =<< ls_convert_to_flat_spin_basis basisPtrPtr normalBasisPtr
      peek basisPtrPtr
  FlatSpinBasis <$> newForeignPtr ls_destroy_flat_spin_basis flatBasisPtr
  where
    (SpinBasis normalBasis) = mkBasis spec

mkBasis :: BasisSpec -> SpinBasis
mkBasis (BasisSpec numberSpins hammingWeight spinInversion symmetries)
  | numberSpins > 64 =
    error $
      "invalid number_spins: " <> show numberSpins <> "; exact diagonalization is not feasible "
        <> "for systems larger than 64 spins"
  | otherwise =
    let symmetryGroup = mkSymmetryGroup (toSymmetry <$> symmetries)
     in System.IO.Unsafe.unsafePerformIO $
          mkRawBasis symmetryGroup numberSpins hammingWeight spinInversion
-}

{-
mkInteraction :: InteractionSpec -> Interaction
mkInteraction (InteractionSpec matrix sites) =
  System.IO.Unsafe.unsafePerformIO $
    mkInteraction' matrix sites

mkOperator :: SpinBasis -> OperatorSpec -> Operator
mkOperator basis (OperatorSpec _ interactions) =
  mkRawOperator basis $
    mkInteraction <$> interactions
-}

{-
mkRawOperator :: SpinBasis -> NonEmpty Interaction -> Operator
mkRawOperator (SpinBasis basis) terms = System.IO.Unsafe.unsafePerformIO $ do
  (code, ptr) <- liftIO $
    alloca $ \ptrPtr -> do
      c <- withInteractions (toList terms) $ \n interactionsPtr ->
        withForeignPtr basis $ \basisPtr ->
          ls_create_operator ptrPtr basisPtr (fromIntegral n) interactionsPtr
      if c == 0
        then (,) <$> pure c <*> peek ptrPtr
        else pure (c, nullPtr)
  checkStatus code
  fmap Operator . liftIO $ newForeignPtr ls_destroy_operator ptr

withInteractions :: [Interaction] -> (Int -> Ptr (Ptr Cinteraction) -> IO a) -> IO a
withInteractions xs func = withManyForeignPtr pointers func
  where
    pointers = (\(Interaction p) -> p) <$> xs
-}

{-
toMatrix ::
  (Monad m, Show r, Real r) =>
  -- | Expected dimension @n@ of the matrix
  Int ->
  -- | Square @matrix@ of dimension @n@
  [[Complex r]] ->
  -- | Row-major representation of @matrix@
  m (Vector (Complex Double))
toMatrix !dim !rows = do
  when (length rows /= dim) . impureThrow . SpinEDException $
    "invalid matrix: " <> show rows <> "; expected a square matrix of dimension " <> show dim
  when (any ((/= dim) . length) rows) . impureThrow . SpinEDException $
    "invalid matrix: " <> show rows <> "; expected a square matrix of dimension " <> show dim
  return . V.fromList . fmap (fmap (fromRational . toRational)) . concat $ rows

type CreateInteraction = Ptr (Ptr Cinteraction) -> Ptr (Complex Double) -> CUInt -> Ptr CUShort -> IO CInt

class MakeInteraction a where
  mkInteraction' :: (MonadIO m, MonadThrow m, Show r, Real r) => [[Complex r]] -> [a] -> m Interaction

instance MakeInteraction Int where
  mkInteraction' matrix sites =
    toMatrix 2 matrix >>= \matrix' ->
      unsafeMkInteraction ls_create_interaction1 (length sites) matrix' sites'
    where
      sites' = V.fromList sites

instance MakeInteraction (Int, Int) where
  mkInteraction' matrix sites =
    toMatrix 4 matrix >>= \matrix' ->
      unsafeMkInteraction ls_create_interaction2 (length sites) matrix' sites'
    where
      sites' = V.fromList $ sites >>= \(x₁, x₂) -> [x₁, x₂]

instance MakeInteraction (Int, Int, Int) where
  mkInteraction' matrix sites =
    toMatrix 8 matrix >>= \matrix' ->
      unsafeMkInteraction ls_create_interaction3 (length sites) matrix' sites'
    where
      sites' = V.fromList $ sites >>= \(x₁, x₂, x₃) -> [x₁, x₂, x₃]

instance MakeInteraction (Int, Int, Int, Int) where
  mkInteraction' matrix sites =
    toMatrix 16 matrix >>= \matrix' ->
      unsafeMkInteraction ls_create_interaction4 (length sites) matrix' sites'
    where
      sites' = V.fromList $ sites >>= \(x₁, x₂, x₃, x₄) -> [x₁, x₂, x₃, x₄]

instance MakeInteraction [Int] where
  mkInteraction' _ [] =
    throwIO . SpinEDException $
      "zero-point interactions (i.e. constant factors) are not supported"
  mkInteraction' matrix rows@(r : _) = case n of
    1 -> mkInteraction' matrix =<< mapM match1 rows
    2 -> mkInteraction' matrix =<< mapM match2 rows
    3 -> mkInteraction' matrix =<< mapM match3 rows
    4 -> mkInteraction' matrix =<< mapM match4 rows
    _ ->
      throwIO . SpinEDException $
        "currently only 1-, 2-, 3-, and 4-point interactions are supported, but received n="
          <> show n
    where
      n = length r
      match1 [x₁] = return x₁
      match1 _ = throwIO failure
      match2 [x₁, x₂] = return (x₁, x₂)
      match2 _ = throwIO failure
      match3 [x₁, x₂, x₃] = return (x₁, x₂, x₃)
      match3 _ = throwIO failure
      match4 [x₁, x₂, x₃, x₄] = return (x₁, x₂, x₃, x₄)
      match4 _ = throwIO failure
      failure =
        SpinEDException $
          "invalid sites: " <> show rows <> "; expected an array of length-" <> show n <> " tuples"
-}

{-
unsafeMkInteraction ::
  (MonadIO m, MonadThrow m) =>
  CreateInteraction ->
  Int ->
  Vector (Complex Double) ->
  Vector Int ->
  m Interaction
unsafeMkInteraction ffiCreate numberSites matrix sites = do
  when (V.any (< 0) sites) . throwIO . SpinEDException $
    "site indices must be all non-negative numbers"
  (code, ptr) <- liftIO $
    alloca $ \ptrPtr -> do
      c <- V.unsafeWith matrix $ \matrixPtr ->
        V.unsafeWith (V.map fromIntegral sites) $ \sitesPtr ->
          ffiCreate ptrPtr matrixPtr (fromIntegral numberSites) sitesPtr
      if c == 0
        then (,) <$> pure c <*> peek ptrPtr
        else pure (c, nullPtr)
  checkStatus code
  fmap Interaction . liftIO $ newForeignPtr ls_destroy_interaction ptr
-}

{-
-- | Extension of 'withForeignPtr' to lists of 'ForeignPtr's.
withManyForeignPtr :: [ForeignPtr a] -> (Int -> Ptr (Ptr a) -> IO b) -> IO b
withManyForeignPtr xs func = loop [] xs
  where
    loop acc (y : ys) = withForeignPtr y $ \y' -> loop (y' : acc) ys
    loop acc [] = withArrayLen (reverse acc) func

withSymmetries :: [Symmetry] -> (Int -> Ptr (Ptr CSymmetry) -> IO a) -> IO a
withSymmetries xs func = withManyForeignPtr pointers func
  where
    pointers = (\(Symmetry p) -> p) <$> xs

mkObject :: (Ptr (Ptr a) -> IO CInt) -> IO (Ptr a)
mkObject f = alloca $ \ptrPtr -> f ptrPtr >>= checkStatus >> peek ptrPtr
-}

  {-
foreign import ccall unsafe "lattice_symmetries/lattice_symmetries.h ls_error_to_string"
  ls_error_to_string :: CInt -> IO CString

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_destroy_string"
  ls_destroy_string :: CString -> IO ()

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_enable_logging"
  ls_enable_logging :: IO ()

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_disable_logging"
  ls_disable_logging :: IO ()

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_create_symmetry"
  ls_create_symmetry :: Ptr (Ptr CSymmetry) -> CUInt -> Ptr CUInt -> CUInt -> IO CInt

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h &ls_destroy_symmetry"
  ls_destroy_symmetry :: FunPtr (Ptr CSymmetry -> IO ())

foreign import ccall unsafe "lattice_symmetries/lattice_symmetries.h ls_create_group"
  ls_create_group :: Ptr (Ptr CGroup) -> CUInt -> Ptr (Ptr CSymmetry) -> IO CInt

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h &ls_destroy_group"
  ls_destroy_group :: FunPtr (Ptr CGroup -> IO ())

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_get_group_size"
  ls_get_group_size :: Ptr CGroup -> IO CUInt

foreign import ccall unsafe "lattice_symmetries/lattice_symmetries.h ls_create_spin_basis"
  ls_create_spin_basis :: Ptr (Ptr Cspin_basis) -> Ptr CGroup -> CUInt -> CInt -> CInt -> IO CInt

foreign import ccall unsafe "lattice_symmetries/lattice_symmetries.h &ls_destroy_spin_basis"
  ls_destroy_spin_basis :: FunPtr (Ptr Cspin_basis -> IO ())

foreign import ccall safe "lattice_symmetries/lattice_symmetries.h ls_build"
  ls_build :: Ptr Cspin_basis -> IO CInt

foreign import ccall safe "lattice_symmetries/lattice_symmetries.h ls_build_unsafe"
  ls_build_unsafe :: Ptr Cspin_basis -> Word64 -> Ptr Word64 -> IO CInt

foreign import ccall unsafe "lattice_symmetries/lattice_symmetries.h ls_get_number_states"
  ls_get_number_states :: Ptr Cspin_basis -> Ptr Word64 -> IO CInt

{- ORMOLU_DISABLE -}
foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_flat_spin_basis_state_info"
  ls_flat_spin_basis_state_info :: Ptr CFlatSpinBasis -> Word64 -> Ptr () -> Ptr () ->
                                   Ptr (Complex Double) -> Ptr Double -> IO ()
{- ORMOLU_ENABLE -}

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_convert_to_flat_spin_basis"
  ls_convert_to_flat_spin_basis :: Ptr (Ptr CFlatSpinBasis) -> Ptr Cspin_basis -> IO CInt

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h &ls_destroy_flat_spin_basis"
  ls_destroy_flat_spin_basis :: FunPtr (Ptr CFlatSpinBasis -> IO ())

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_create_interaction1"
  ls_create_interaction1 :: Ptr (Ptr Cinteraction) -> Ptr (Complex Double) -> CUInt -> Ptr CUShort -> IO CInt

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_create_interaction2"
  ls_create_interaction2 :: Ptr (Ptr Cinteraction) -> Ptr (Complex Double) -> CUInt -> Ptr CUShort -> IO CInt

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_create_interaction3"
  ls_create_interaction3 :: Ptr (Ptr Cinteraction) -> Ptr (Complex Double) -> CUInt -> Ptr CUShort -> IO CInt

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_create_interaction4"
  ls_create_interaction4 :: Ptr (Ptr Cinteraction) -> Ptr (Complex Double) -> CUInt -> Ptr CUShort -> IO CInt

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h ls_interaction_is_real"
  ls_interaction_is_real :: Ptr () -> IO CBool

foreign import capi unsafe "lattice_symmetries/lattice_symmetries.h &ls_destroy_interaction"
  ls_destroy_interaction :: FunPtr (Ptr Cinteraction -> IO ())

foreign import ccall unsafe "lattice_symmetries/lattice_symmetries.h ls_create_operator"
  ls_create_operator :: Ptr (Ptr Coperator) -> Ptr Cspin_basis -> CUInt -> Ptr (Ptr Cinteraction) -> IO CInt

foreign import ccall unsafe "lattice_symmetries/lattice_symmetries.h &ls_destroy_operator"
  ls_destroy_operator :: FunPtr (Ptr Coperator -> IO ())
-}

max_stars_repo_name: twesterhout/lattice-symmetries-haskell

max_issues_repo_name: twesterhout/lattice-symmetries-haskell

max_forks_repo_name: twesterhout/lattice-symmetries-haskell


text:
module Math.FROG.Tools
( mkTrace
, mkSigGate
, shift
, mkInitialGuess
, flatPolar2Complex
, center
) where

import qualified System.Random as R
import qualified Numeric.LinearAlgebra as LA
import qualified Numeric.GSL.Fourier as F
import qualified Foreign

import Math.FROG.Types

mkTrace :: ComplexSignal -> ComplexSignal -> Trace
mkTrace field gate = LA.fromColumns $ map atEachDelayApply [upper,upper-1..lower]
    where
    atEachDelayApply = (^2) . (LA.mapVector LA.magnitude)
                            . (flip shift upper)
                            .  F.fft
                            . (*field)
                            . (shift gate)
    sz = LA.dim field
    upper = (sz `quot` 2)
    lower = (-1) * (sz `quot` 2) + 1

mkSigGate :: Nonlinearity -> ComplexSignal -> ComplexSignal -> (ComplexSignal, ComplexSignal)
mkSigGate PG f g = (id f, LA.mapVector ((^2) . abs) g)
mkSigGate SHG f g = (id f, id g)
mkSigGate SD f g = (LA.mapVector (^2) f, LA.mapVector LA.conjugate g)
mkSigGate THG f g = (LA.mapVector (^2) f, id g)

shift :: ComplexSignal -> Int -> ComplexSignal
shift gate delay = LA.buildVector sz (builder delay)
    where
    sz = LA.dim gate
    builder d k = if any (<0) [(sz-1-d-k), (k+d)] then (getElem (k+d-(signum d)*sz)) else (getElem (k+d))
        where
        getElem = (LA.@>) gate

mkInitialGuess :: Int -> IO (LA.Vector Double)
mkInitialGuess sz = do
    stdgen <- R.getStdGen
    let randMags = R.randomRs (0, 1.0) stdgen
    let randAngs = R.randomRs (-pi, pi) stdgen
    return $ LA.buildVector (2*sz) (\k -> if k < sz then (randMags !! k) else (randAngs !! k))

flatPolar2Complex :: LA.Vector Double -> LA.Vector (LA.Complex Double)
flatPolar2Complex v = LA.buildVector sz (\k -> LA.mkPolar (mags LA.@> k) (angs LA.@> k))
    where
    sz = (LA.dim v) `quot` 2
    mags = LA.subVector 0 sz v
    angs = LA.subVector sz sz v

center :: ComplexSignal -> ComplexSignal
center f = shift f $ getDelay . getMaxInd $ f
    where
    sz = LA.dim f
    getMaxInd = LA.maxIndex . (LA.mapVector LA.magnitude)
    getDelay = (*(-1)) . ((-) (sz `quot` 2))


text:
module Scheme.Eval.Primitives (ioPrimitives, basicPrimitives) where

import Control.Monad
import Data.Maybe
import Data.Complex as C
import Data.Ratio
import Control.Monad.Error
import Scheme.Types
import System.IO

import Scheme.Eval.DynamicTypes

prim = undefined

convert typ wrap = expect typ >=> return . wrap

foldlOnType typ f = foldM (\acc e -> if check typ e
    then return (acc `f` e)
    else errorNeedType typ (show e))
foldlOnType1 typ f initPrep l = do expect typ (head l)
                                   foldlOnType typ f (initPrep $ head l) (init l)

basicPrimitives :: [(String, [LispVal] -> ThrowsErrorIO LispVal)]
basicPrimitives =
    let toBool = return . Bool
        dol = dottedOrPlainListType
    in
    -- number type checking
    [ ("number?",
        oneArgOnly >=> toBool . isNum)
    , ("real?",
        oneArgOnly >=> toBool . checkNum isDouble)
    , ("rational?",
        oneArgOnly >=> toBool . checkNum isRational)
    , ("integer?",
        oneArgOnly >=> toBool . checkNum isInt)
    , ("complex?",
        oneArgOnly >=> toBool . checkNum isComplex)

    -- only nice rationals and integers are exact
    , ("exact?",
        oneArgOnly >=> \n -> toBool $ any ($n) [checkNum isRational, checkNum isInt])

    -- arithmetic, preserve exactness in this order: integer, rational, float, complex
    , ("+",
        minArgs 2 >=> expect (listOf numType) . List >=> return . Num . sum)
    , ("-",
        minArgs 2 >=> expect (listOf numType) . List >=> return . Num . foldl1 (-))
    , ("*",
        minArgs 2 >=> expect (listOf numType) . List >=> return . Num . product)
    , ("/",
        onlyArgs 2 >=> \[_a, _b] -> do a <- expect numType _a
                                       b <- expect numType _b
                                       return . Num $ a / b)
    , ("negate",
        oneArgOnly >=> expect numType >=> return . Num . negate)
    , ("signum",
        oneArgOnly >=> expect numType >=> return . Num . signum)
    , ("abs",
        oneArgOnly >=> expect numType >=> return . Num . abs)
    , ("recip",
        oneArgOnly >=> expect numType >=> return . Num . recip)

    -- integer arithmetic functions
    , ("mod",
        minArgs 2 >=> foldlOnType1 numIntType (\(Num (Int a)) (Num (Int b))
                                              -> Num . Int $ a `mod` b) id)
    , ("quotient",
        minArgs 2 >=> foldlOnType1 numIntType (\(Num (Int a)) (Num (Int b))
                                              -> Num . Int $ a `quot` b) id)
    , ("remainder",
        minArgs 2 >=> foldlOnType1 numIntType (\(Num (Int a)) (Num (Int b))
                                              -> Num . Int $ a `rem` b) id)

    -- boolean
    , ("boolean?", oneArgOnly >=> toBool . isBool)
    , ("&&",
        minArgs 1 >=> liftM Bool . foldlOnType boolType (\a b -> a && fromJust (getBool b)) True)
    , ("||",
        minArgs 1 >=> liftM Bool . foldlOnType boolType (\a b -> a || fromJust (getBool b)) False)

    -- list processing !
    , ("list?",
        oneArgOnly >=> \l -> toBool $ any ($l) [isList, isDottedList, isVector])
    , ("dottedlist?",
        oneArgOnly >=> toBool . isDottedList)

    , ("car",
        oneArgOnly >=> expect dol >=> return .
                        (return . either id fst >=> safeHead)
                   >=> maybe (throwError $ Default "head-less list given to car") return)

    , ("cdr",
        oneArgOnly >=> expect dol >=> return .
                        --extract the dotted/normal list
                        either (liftM List . safeTail)
                            -- handle two cases for dotted lists
                            (\(xs, dot) -> Just $ if null xs then dot else DottedList (tail xs) dot)

                    -- handle bad tail
                   >=> maybe (throwError $ Default "tail-less list given to cdr") return)

    , ("cons", onlyArgs 2 >=> (\(a:b:[]) ->
                maybe (return $ List [a,b]) return
                    (get dol b >>= either (Just . List . (a:)) (\(xs, dot) -> Just $ DottedList (a:xs) dot))
                ))

    , ("vector?",
        oneArgOnly >=> toBool . isVector)
    , ("vector-ref", prim)

    , ("vector-set", prim)
    , ("list->vector",
        oneArgOnly >=> prim)
    , ("vector->list",
        oneArgOnly >=> prim)

    -- symbol processing
    , ("symbol?",
        oneArgOnly >=> toBool . isAtom)

    , ("string->symbol",
        oneArgOnly >=> convert stringType Atom)
    , ("symbol->string",
        oneArgOnly >=> convert atomType String)

    -- string processing
    , ("string?",
        oneArgOnly >=> toBool . isString)
    , ("string-ref",
        onlyArgs 2 >=> \([lispstr, lispix]) -> do
            str <- expect stringType lispstr
            ix  <- expect numIntType lispix
            char           <- maybe (throwError $ BadExpr "index out of bounds" lispix) return (safeIndex (fromIntegral ix) str)
            return . Character $ char)

    -- string construction
    , ("string->list",
        oneArgOnly >=> convert stringType (List . map Character))
    , ("list->string",
        convert (listOf charType) String . List)

    -- characters
    , ("character?",
        oneArgOnly >=> toBool . isCharacter)
    , ("character=", minArgs 1 >=> liftM (Bool . fst) . foldlOnType charType (\(res, prev) c -> (isNothing prev || getCharacter c == prev, getCharacter c)) (True, Nothing))

    -- equivalence
    , ("eq?",    onlyArgs 2 >=> return . eqv)
    , ("eqv?",   onlyArgs 2 >=> return . eqv)
    , ("=",     onlyArgs 2 >=> return . eqv)
    , ("equal?", onlyArgs 2 >=> return . equal)
    , ("==",    onlyArgs 2 >=> return . equal)
    ]

ioPrimitives =
    [ ("open-input-file", oneArgOnly >=> expect stringType >=> makePort ReadMode)
    , ("open-output-file", oneArgOnly >=> expect stringType >=> makePort WriteMode)
    , ("close-port", closePort)
    , ("close-input-port", closePort)
    , ("close-output-port", closePort)
    , ("read-contents", oneArgOnly >=> expect portType >=> liftM String . liftIO . hGetContents)
    , ("write", maxArgs 2 >=> minArgs 1 >=> \l -> case l of
        [obj] -> safeLiftIO (print obj) >> lispNull
        [obj, p] -> expect portType p >>= \port -> safeLiftIO (hPrint port (show obj)) >> lispNull)
    ]

closePort = oneArgOnly >=> liftIO . hClose . fromJust . getPort >=> const lispNull

makePort mode = liftM Port . safeLiftIO . flip openFile mode

eqv [a,b]
    | a == List [] && a == b = Bool True
    | all (\x -> check (dottedListType `orType` listType `orType` vectorType) x) [a,b] = Bool False
    | otherwise = equal [a,b]

equal [a,b] = Bool $ a == b

max_stars_repo_name: oroguh/SCMinHS

max_issues_repo_name: oroguh/SCMinHS

max_forks_repo_name: oroguh/SCMinHS


text:
module MachineLearning.MultiSvmClassifierTest
(
  tests
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Approx
import Test.HUnit.Plus

import MachineLearning.DataSets (dataset2)

import qualified Numeric.LinearAlgebra as LA
import qualified MachineLearning as ML
import MachineLearning.Optimization
import MachineLearning.Model
import MachineLearning.MultiSvmClassifier

(x, y) = ML.splitToXY dataset2

model = MultiClass (MultiSvm 1 2)

x1 = ML.addBiasDimension x
onesTheta :: LA.Vector LA.R
onesTheta = LA.konst 1 (2 * LA.cols x1)
zeroTheta :: LA.Vector LA.R
zeroTheta = LA.konst 0 (2 * LA.cols x1)

processX muSigma x = ML.addBiasDimension $ ML.featureNormalization muSigma $ ML.mapFeatures 6 x

muSigma = ML.meanStddev (ML.mapFeatures 6 x)
x2 = processX muSigma x

xPredict = LA.matrix 2 [ -0.5, 0.5
                       , 0.2, -0.2
                       , 1, 1
                       , 1, 0
                       , 0, 0
                       , 0, 1]
xPredict2 = processX muSigma xPredict
yExpected = LA.vector [1, 1, 0, 0, 1, 0]

gradientCheckingEps :: Double
gradientCheckingEps = 1e-3

eps = 0.0001

zeroTheta2 = LA.konst 0 (2 * LA.cols x2)
(thetaGD, _) = minimize (GradientDescent 0.001) model eps 150 (L2 1) x2 y zeroTheta2
(thetaCGFR, _) = minimize (ConjugateGradientFR 0.1 0.1) model eps 30 (L2 0.5) x2 y zeroTheta2
(thetaCGPR, _) = minimize (ConjugateGradientPR 0.1 0.1) model eps 30 (L2 0.5) x2 y zeroTheta2
(thetaBFGS, _) = minimize (BFGS2 0.1 0.1) model eps 30 (L2 0.5) x2 y zeroTheta2

checkGradientTest lambda theta eps = do
  let diffs = take 5 $ map (\e -> checkGradient model lambda x1 y theta e) [1e-3, 1.1e-3 ..]
      diff = minimum diffs
  assertApproxEqual "" eps 0 diff

tests = [  testGroup "gradient checking" [
            testCase "non-zero theta, non-zero lambda" $ checkGradientTest (L2 2) onesTheta 3e-2
              , testCase "zero theta, non-zero lambda" $ checkGradientTest (L2 2) zeroTheta gradientCheckingEps
              , testCase "non-zero theta, zero lambda" $ checkGradientTest (L2 0) onesTheta gradientCheckingEps
              , testCase "zero theta, zero lambda" $ checkGradientTest (L2 0) zeroTheta gradientCheckingEps
              , testCase "non-zero theta, no reg" $ checkGradientTest RegNone onesTheta gradientCheckingEps
              , testCase "zero theta, no reg" $ checkGradientTest RegNone zeroTheta gradientCheckingEps
              ]
        , testGroup "learn" [
            testCase "Gradient Descent" $ assertVector "" 0.01 yExpected (hypothesis model xPredict2 thetaGD)
            , testCase "Conjugate Gradient FR" $ assertVector "" 0.01 yExpected (hypothesis model xPredict2 thetaCGFR)
            , testCase "Conjugate Gradient PR" $ assertVector "" 0.01 yExpected (hypothesis model xPredict2 thetaCGPR)
            , testCase "BFGS" $ assertVector "" 0.01 yExpected (hypothesis model xPredict2 thetaBFGS)
            ]
        ]

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Main where
import           Approx
import           Data.Complex
import           Carlson
import           Elliptic
import           Test.Tasty       (defaultMain, testGroup)
import           Test.Tasty.HUnit (assertEqual, testCase)

i :: Complex Double
i = 0.0 :+ 1.0

main :: IO ()
main = defaultMain $
  testGroup "Tests"
  [ testCase "RF value 1" $
      assertEqual ""
        (approx 12 (carlsonRF 1 2 0))
        (approx 12 1.3110287771461),

    testCase "RF value 2" $
      assertEqual ""
        (approx 12 (carlsonRF i (-i) 0))
        (approx 12 1.8540746773014),

    testCase "RF value 3" $
      assertEqual ""
        (approx 12 (carlsonRF 0.5 1 0))
        (approx 12 1.8540746773014),

    testCase "RF value 4" $
      assertEqual ""
        (approx 13 (carlsonRF (i-1) i 0))
        (approx 13 (0.79612586584234 :+ (-1.2138566698365))),

    testCase "RF value 5" $
      assertEqual ""
        (approx 13 (carlsonRF 2 3 4))
        (approx 13 0.58408284167715),

    testCase "RF value 6" $
      assertEqual ""
        (approx 12 (carlsonRF i (-i) 2))
        (approx 12 1.0441445654064),

    testCase "RF value 7" $
      assertEqual ""
        (approx 13 (carlsonRF (i-1) i (1-i)))
        (approx 13 (0.93912050218619 :+ (-0.53296252018635))),

    testCase "RC value 1" $
      assertEqual ""
        (approx 14 (carlsonRC 0 0.25))
        (approx 14 pi),

    testCase "RC value 2" $
      assertEqual ""
        (approx 14 (carlsonRC 2.25 2))
        (approx 14 (log 2)),

    testCase "RC value 3" $
      assertEqual ""
        (approx 12 (carlsonRC 0 i))
        (approx 12 ((1-i)*1.1107207345396)),

    testCase "RC value 4" $
      assertEqual ""
        (approx 13 (carlsonRC (-i) i))
        (approx 13 (1.2260849569072 :+ (-0.34471136988768))),

    testCase "RC value 5" $
      assertEqual ""
        (approx 14 (carlsonRC 0.25 (-2)))
        (approx 14 ((log 2 :+ (-pi))/ 3)),

    testCase "RC x y = RF x y y" $ do
      let x = 5 :+ 6
          y = 2 :+ (-9)
      assertEqual ""
        (approx 14 (carlsonRC x y))
        (approx 14 (carlsonRF x y y)),

    testCase "RJ x y y p" $ do
      let x = 1 :+ 1
          y = (-2) :+ 3
          p = 0 :+ 4
      assertEqual ""
        (approx 14 (carlsonRJ x y y p))
        (approx 14 (3*(carlsonRC x y - carlsonRC x p) / (p-y))),

    testCase "RJ homogeneity" $ do
      let x = 1 :+ 1
          y = (-2) :+ 3
          z = -3
          p = 0 :+ 4
          kappa = 2 :+ 0
      assertEqual ""
        (approx 14 (carlsonRJ x y z p / kappa / sqrt kappa))
        (approx 14 (carlsonRJ (kappa*x) (kappa*y) (kappa*z) (kappa*p))),

    testCase "Complete elliptic integral K" $ do
      let m = 2 :+ (-3)
      assertEqual ""
        (approx 14 (ellipticF (pi/2) m))
        (approx 14 (carlsonRF 0 (1-m) 1)),

    testCase "Complete ellipticE - RG" $ do
      let m = 2 :+ (-3)
      assertEqual ""
        (approx 14 (ellipticE (pi/2) m))
        (approx 14 (2 * carlsonRG 0 (1-m) 1)),

    testCase "Complete ellipticE - RD" $ do
      let m = 2 :+ (-3)
      assertEqual ""
        (approx 14 (ellipticE (pi/2) m))
        (approx 14 ((1-m) * (carlsonRD 0 (1-m) 1 + carlsonRD 0 1 (1-m)) / 3)),

    testCase "jacobiZeta m=1" $ do
      let z = (-1) :+ 8
      assertEqual ""
        (approx 14 (jacobiZeta z 1))
        (approx 14 (sin z))

  ]

max_stars_repo_name: stla/elliptic

max_issues_repo_name: stla/elliptic

max_forks_repo_name: stla/elliptic


text:
import Control.Arrow
import System.Directory
import Data.Char
import qualified Data.Text as T
import qualified Data.Text.IO as T
import Data.List
import qualified Data.Set as S
import qualified Data.Map.Strict as M
import qualified Data.IntMap.Strict as IM
import NLP.Stemmer
import Data.Ord
import Data.Function
import Data.SVM
import qualified Numeric.LinearAlgebra.Data as LA
import qualified Numeric.LinearAlgebra as LA

-- Ratings are embedded in the file names, e.g., 100_7.txt has a rating of 7/10. Extract it.
getRating :: (Fractional a, Read a) => String -> a
getRating = (\x -> (x-5.5)/4.5) . read . takeWhile (/= '.') . drop 1 . dropWhile (/= '_')

-- Get the rating and text from all files in a directory.
readSamples :: (Fractional a, Read a) => String -> IO [(a, T.Text)]
readSamples dir = mapM ((sequence .) $ getRating &&& T.readFile . (dir++)) . filter ((/= '.') . head) =<< getDirectoryContents dir

-- All of the training data, which is located across two directories.
trainingSamples :: (Fractional a, Read a) => IO [(a, T.Text)]
trainingSamples = do
  pos <- readSamples "aclImdb/train/pos/"
  neg <- readSamples "aclImdb/train/neg/"
  return $ takeEvery 5 (pos ++ neg)
validationSamples :: (Fractional a, Read a) => IO [(a, T.Text)]
validationSamples = do
  pos <- readSamples "aclImdb/train/pos/"
  neg <- readSamples "aclImdb/train/neg/"
  return $ takeEvery 10 (undefined:(pos ++ neg))

-- Words to ignore
stoplist :: S.Set T.Text
stoplist = S.fromList $ map T.pack ["i", "me", "my", "she", "her", "he", "him", "his", "they", "them", "their", "it", "who", "this", "that", "which",
  "", "br", "the", "a", "an", "is", "was", "be", "have", "are", "had", "were", "do", "did",
  "and", "or", "to", "in", "for", "with", "as", "one", "at", "by", "about", "from", "other", "into", "when", "would", "than",
  "also", "how", "becaus", "could", "after", "thing", "dont", "year", "two", "after", "peopl", "get", "where", "did", "off"]

minimalStoplist = S.fromList $ map T.pack ["the", "a", "an", "to", "is", "that", ""]
negateList = S.fromList $ map T.pack ["not", "isnt"]

-- Word splitter which ignores puctuation and doesn't care about capitalization or apostrophes.
words' :: T.Text -> [T.Text]
words' = map (stem' English) . filter (not . T.null) . T.split (not . isAlpha) . T.toLower . T.filter (/= '\'')

data PreFilterMode = Flip | Drop | Pass

-- Apply +/- modifier to word list
preFilter :: Num a => PreFilterMode -> [T.Text] -> [(a, T.Text)]
preFilter _ [] = []
preFilter Flip (w:ws)
  | w `S.member` negateList = alterHead (first negate) (preFilter Flip ws)
  | w `S.member` minimalStoplist = preFilter Flip ws
  | T.pack "un" `T.isPrefixOf` w = (-1, T.drop 2 w) : preFilter Flip ws
  | T.pack "nt" `T.isSuffixOf` w = (-1, T.dropEnd 2 w) : preFilter Flip ws
  | otherwise = (1, w) : preFilter Flip ws
  where alterHead f [] = []
        alterHead f (x:xs) = (f x) : xs
preFilter Drop (w:ws)
  | w `S.member` negateList = preFilter Drop ws
  | w `S.member` minimalStoplist = preFilter Drop ws
  | T.pack "un" `T.isPrefixOf` w = (1, T.drop 2 w) : preFilter Drop ws
  | T.pack "nt" `T.isSuffixOf` w = (1, T.dropEnd 2 w) : preFilter Drop ws
  | otherwise = (1, w) : preFilter Drop ws
preFilter Pass ws = zip (repeat 1) ws

-- Unify the number of occurences
bag :: Num a => [(a, T.Text)] -> [(a, T.Text)]
bag = map (sum . map fst &&& snd . head) . groupBy (curry $ uncurry (==) . both snd) . sortBy (comparing snd)

-- The stemming method from NLP.Stemmer, wrapped up to use Text insead of String
stem' :: Stemmer -> T.Text -> T.Text
stem' lang = T.pack . stem lang . T.unpack

-- The set of words in a dataset that meet the common threshold
common :: PreFilterMode -> Int -> [(a, T.Text)] -> S.Set T.Text
common mode thresh dataset = S.fromList . map snd . filter (\w -> fst w * thresh >= datalen) $ (df mode (`S.notMember` stoplist) dataset)
  where datalen = length dataset

-- Number of documents with each word not excluded by filter
df :: PreFilterMode -> (T.Text -> Bool) -> [(a, T.Text)] -> [(Int, T.Text)]
df mode wf = map (length &&& head) . group . sort . concatMap (map head . group . sort . filter wf . map snd . preFilter mode . words' . snd)

df' :: PreFilterMode -> (T.Text -> Bool) -> [(a, T.Text)] -> [(Float, T.Text)]
df' mode wf = map (log . balance . sum . map fst &&& snd . head) .
              groupBy ((==) `on` snd) . sortBy (comparing snd) .
              concatMap (bag . filter (wf . snd) . preFilter mode . words' . snd)

-- Calculates the Delta Inverse Document Frequency of the dataset
didf :: (Num a, Ord a) => PreFilterMode -> S.Set T.Text -> [(a, T.Text)] -> M.Map T.Text Float
didf mode legalWords dataset = uncurry (M.unionWith (+)) . both (M.fromList . map swap) .
    (                     df' mode (`S.member` legalWords) . filter ((>0) . fst) &&&
     map (first negate) . df' mode (`S.member` legalWords) . filter ((<0) . fst)) $ dataset

-- Converts a review into a vector for linear algebra
t2v :: PreFilterMode -> M.Map T.Text Int -> M.Map T.Text Double -> T.Text -> LA.Vector Double
t2v mode lut weights t = LA.fromList $ map (flip (IM.findWithDefault 0) (t2v' mode lut weights t)) [0 .. length lut - 1]

-- Converts a review into a vector for SVM, which is a different format
t2v' :: PreFilterMode -> M.Map T.Text Int -> M.Map T.Text Double -> T.Text -> IM.IntMap Double
t2v' mode lut weights t = IM.fromList . map (\(x, w) -> (lut M.! w, realToFrac $ weights M.! w * x)) .
                          filter ((`M.member` lut) . snd) . bag . preFilter mode . words' $ t

-- Converts a dataset into a term frequency matrix for linear algebra
data2m :: PreFilterMode -> M.Map T.Text Int -> M.Map T.Text Double -> [(a, T.Text)] -> LA.Matrix Double
data2m mode lut weights = LA.fromRows . map (t2v mode lut weights . snd)

-- Converts a matrix into a SVM problem
m2p :: [(Double, a)] -> LA.Matrix Double -> [(Double, IM.IntMap Double)]
m2p dataset = zip (map (clamp . fst) dataset) . map (\row -> IM.fromList $ map (id &&& (row LA.!)) [0 .. LA.size row - 1]) . LA.toRows

-- Applies a SVM model to a problem
checkProblem :: Model -> [(Double, IM.IntMap Double)] -> IO [(Double, Double)]
checkProblem model = sequence . map (sequence . second (predict model))

showReport :: (Ord a, Ord b, Num a, Num b) => [(a, b)] -> String
showReport results = show a ++ "/" ++ show n ++ " (" ++ show (100 * fromIntegral a / fromIntegral n) ++ "%)"
  where a = (length $ filter (\(a, b) -> (a>0) == (b>0)) results)
        n = (length results)

printReport :: (Ord a, Ord b, Num a, Num b) => (c, [(a, b)], [(a, b)]) -> IO ()
printReport (_, trainingResults, testingResults) = do
  putStrLn $ "\ttraining accuracy: " ++ showReport trainingResults ++ "\tvalidation accuracy: " ++ showReport testingResults

runExperiment :: (PreFilterMode, Int, Int, KernelType, Algorithm, [(Double, T.Text)], [(Double, T.Text)]) ->
                 IO (S.Set T.Text, [(Double, Double)], [(Double, Double)])
runExperiment (filterMode, commonThreshold, nSingularValues, svmKernel, svmAlgorithm, trainingDataset, testingDataset) = do
  let legalWords = common filterMode commonThreshold trainingDataset
  let weights = realToFrac <$> didf filterMode legalWords trainingDataset
  let wordIndexes = M.fromList $ zip (S.toList legalWords) [0..]
  let (u, s, v) = LA.thinSVD (data2m filterMode wordIndexes weights trainingDataset)
  let u' = u LA.?? (LA.All, LA.Take nSingularValues)
  let s' = LA.subVector 0 nSingularValues s
  let v' = v LA.?? (LA.All, LA.Take nSingularValues)
  let testingDataset' = m2p testingDataset $ data2m filterMode wordIndexes weights testingDataset LA.<> v'
  let trainingDataset' = m2p trainingDataset (u' LA.<> LA.diag s')
  svmModel <- withPrintFn (const $ return ()) $ train svmAlgorithm svmKernel trainingDataset'
  resultsTraining <- checkProblem svmModel trainingDataset'
  resultsValidation <- checkProblem svmModel testingDataset'
  return (legalWords, resultsTraining, resultsValidation)

main :: IO ()
main = do
  let thresh = 200
  let nsv = 300
  let kern = Linear --RBF 0.01
  let alg = CSvc 1

  dataset <- trainingSamples
  validationDataset <- validationSamples

  putStrLn "Our Method (Term Frequency Negation):"
  printReport =<< runExperiment (Flip, thresh, nsv, kern, alg, dataset, validationDataset)
  putStrLn "Stoplist & Stem Modifiers Without Negation:"
  printReport =<< runExperiment (Drop, thresh, nsv, kern, alg, dataset, validationDataset)
  putStrLn "Modifiers Treated Normally:"
  printReport =<< runExperiment (Pass, thresh, nsv, kern, alg, dataset, validationDataset)

-- Assorted utility functions

crop :: (Num a, Ord a) => a -> a
crop = max (-1) . min 1

clamp :: (Num a, Ord a) => a -> a
clamp x = if x > 0 then 1 else -1

balance :: (Ord a, Fractional a) => a -> a
balance x = if x > 0 then 1 + x else 1 / (1 - x)

takeEvery :: Int -> [a] -> [a]
takeEvery n xs = case drop (n-1) xs of (y:ys) -> y : takeEvery n ys; [] -> []

dropEvery :: Int -> [a] -> [a]
dropEvery n xs = case splitAt (n-1) xs of (x, (y:ys)) -> x ++ dropEvery n ys; (x, []) -> x

r2 :: Floating a => [(a, a)] -> (a, (a, a))
r2 points = let n = fromIntegral (length points)
                ex = sum (map fst points) / n
                ey = sum (map snd points) / n
                sx = sqrt $ sum (map ((^2) . ((-)ex) . fst) points) / n
                sy = sqrt $ sum (map ((^2) . ((-)ey) . snd) points) / n
                cov = sum (map (\(x, y) -> (x - ex) * (y - ey)) points) / n
            in (cov / (sx * sy), (sx / sy, ex * sy / sx - ey))

both :: (a -> b) -> (a, a) -> (b, b)
both f (x, y) = (f x, f y)

swap :: (a, b) -> (b, a)
swap (x, y) = (y, x)

max_stars_repo_name: LyraSolomon/sentiment-analysis

max_issues_repo_name: LyraSolomon/sentiment-analysis

max_forks_repo_name: LyraSolomon/sentiment-analysis


text:
{-# OPTIONS_GHC -Wall #-}

{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE OverloadedStrings #-}

module Data.NutMeg
    ( NutMeg (..)
    , NutPlot (..)
    , NutWave (..)
    , NutPlotType (..)
    , NutField (..)
    , nutFieldName, nutFieldName'
    , bytesPerReal, bytesPerReal'
    , bytesPerComplex, bytesPerComplex'
    , popNutElement
    , readNutRealRow, readNutComplexRow
    , nutPlot, nutFlag, nutWave
    , nutRealWave, asRealVector
    , nutComplexWave, asComplexVector
    , flattenRealPlots
    , flattenComplexPlots
    , parseNutPlot
    , parseNutMeg
    , readNutRaw
    -- , encodeNutPlot
    ) where

import           Data.Maybe
import           Data.Int
import           Data.Complex
import           Data.Binary.Get
import qualified Data.Map.Strict       as M
import qualified Data.ByteString       as BS
import qualified Data.ByteString.Char8 as CS
import qualified Data.ByteString.Lazy  as BL
import qualified Data.Vector.Unboxed   as V
import qualified Data.Matrix.Unboxed   as A

-- | NutMeg wave form data dypes
data NutWave = NutRealWave    {-# UNPACK #-} !(V.Vector Double)
             | NutComplexWave {-# UNPACK #-} !(V.Vector (Complex Double))
    deriving (Show, Eq)

nutWave :: NutWave -> Either (V.Vector (Complex Double)) (V.Vector Double)
nutWave (NutRealWave nw) = Right nw
nutWave (NutComplexWave nw) = Left nw

nutRealWave :: NutWave -> Maybe (V.Vector Double)
nutRealWave (NutRealWave nw) = Just nw
nutRealWave _                = Nothing

asRealVector :: NutWave -> V.Vector Double
asRealVector (NutRealWave nw) = nw
asRealVector _                = undefined

concatRealWaves :: [NutWave] -> NutWave
concatRealWaves = NutRealWave . V.concat . map asRealVector

joinRealWaves :: NutWave -> NutWave -> NutWave
joinRealWaves a b = concatRealWaves [a, b]

flattenRealPlots :: [NutPlot] -> NutPlot
flattenRealPlots nps = let { nutPlotName  = pn'
                           ; nutNumVars   = nv'
                           ; nutNumPoints = np'
                           ; nutVariables = vn'
                           ; nutPlotType  = pt'
                           ; nutData      = dt'
                           } in NutPlot { .. }
  where fp = head nps
        pn' = nutPlotName fp
        nv' = nutNumVars fp
        np' = nutNumPoints fp
        vn' = nutVariables fp
        pt' = nutPlotType fp
        dt' = M.unionsWith joinRealWaves (map nutData nps)

nutComplexWave :: NutWave -> Maybe (V.Vector (Complex Double))
nutComplexWave (NutComplexWave nw) = Just nw
nutComplexWave _                   = Nothing

asComplexVector :: NutWave -> V.Vector (Complex Double)
asComplexVector (NutComplexWave nw) = nw
asComplexVector _                   = undefined

concatComplexWaves :: [NutWave] -> NutWave
concatComplexWaves = NutComplexWave . V.concat . map asComplexVector

joinComplexWaves :: NutWave -> NutWave -> NutWave
joinComplexWaves a b = concatComplexWaves [a, b]

flattenComplexPlots :: [NutPlot] -> NutPlot
flattenComplexPlots nps = let { nutPlotName  = pn'
                              ; nutNumVars   = nv'
                              ; nutNumPoints = np'
                              ; nutVariables = vn'
                              ; nutPlotType  = pt'
                              ; nutData      = dt'
                              } in NutPlot { .. }
  where fp = head nps
        pn' = nutPlotName fp
        nv' = nutNumVars fp
        np' = nutNumPoints fp
        vn' = nutVariables fp
        pt' = nutPlotType fp
        dt' = M.unionsWith joinRealWaves (map nutData nps)

-- | Real and Complex NutMeg plots
data NutPlotType = NutRealPlot | NutComplexPlot
    deriving (Show, Eq)

-- | Convert a flag from file to plot type
nutFlag :: String -> NutPlotType
nutFlag "complex" = NutComplexPlot
nutFlag "real"    = NutRealPlot
nutFlag _         = undefined

-- | Data Fields in NutMeg Files
data NutField = NutTitle
              | NutDate
              | NutPlotname
              | NutFlags
              | NutNoVariables
              | NutNoPoints
              | NutVariables
              | NutBinData
    deriving (Show, Eq, Enum, Ord, Bounded)

-- Alternaitve Show instance calling `nutFieldName`.
-- instance Show NutField where
--   show nfn = CS.unpack . nutFieldName $ nfn

-- | Translate fieldname to ByteString representation in NutMeg file.
nutFieldName :: NutField -> BS.ByteString
nutFieldName NutTitle       = "Title:"
nutFieldName NutDate        = "Date:"
nutFieldName NutPlotname    = "Plotname:"
nutFieldName NutFlags       = "Flags:"
nutFieldName NutNoVariables = "No. Variables:"
nutFieldName NutNoPoints    = "No. Points:"
nutFieldName NutVariables   = "Variables:"
nutFieldName NutBinData     = "Binary:"

-- | Convenient Transformer to Lazy ByteString
nutFieldName' :: NutField -> BL.ByteString
nutFieldName' = BL.fromStrict . nutFieldName

-- | A 'plot' as it is represented within the NutMeg File format.
data NutPlot = NutPlot { nutPlotName  :: String
                       , nutNumVars   :: Int
                       , nutNumPoints :: Int
                       , nutVariables :: [String]
                       , nutPlotType  :: NutPlotType
                       , nutData      :: !(M.Map String NutWave) }
--              | ComplexPlot { nutPlotName    :: String
--                            , nutNumVars     :: Int
--                            , nutNumPoints   :: Int
--                            , nutVariables   :: [String]
--                            , nutPlotType    :: NutPlotType
--                            , complexNutData :: !(M.Map String NutWave) }
    deriving (Show)

-- | Representation of NutMeg file contents, where plot names are mapped to
-- corresponding plot types and data.
data NutMeg = NutMeg { nutTitle :: String
                     , nutDate  :: String
                     , nutPlots :: !(M.Map String NutPlot)}
    deriving (Show)

-- | How many bytes per real data point
bytesPerReal :: Int
bytesPerReal = 8
-- | Convenient transformation to Int64 for certain calculations
bytesPerReal' :: Int64
bytesPerReal' = fromIntegral bytesPerReal

-- | How many bytes per complex data point
bytesPerComplex :: Int
bytesPerComplex = 2 * bytesPerReal
-- | Convenient transformation to Int64 for certain calculations
bytesPerComplex' :: Int64
bytesPerComplex' = fromIntegral bytesPerReal

-- | Convenience function for reading a NutMeg File
readNutRaw :: FilePath -> IO BS.ByteString
readNutRaw = BS.readFile

-- | Reads a line from a list of ByteStrings, that starts with a given prefix.
readNutElement' :: BS.ByteString -> [BS.ByteString] -> Maybe String
readNutElement' eid bs = CS.unpack . CS.strip <$> elem' bs
    where elem' = CS.stripPrefix eid . head . filter (BS.isPrefixOf eid)

-- | Takes a `NutField` name and returns the correspinding value.
readNutElement :: NutField -> BS.ByteString -> String
readNutElement nf bs = fromJust $ readNutElement' nfn bls
    where nfn = nutFieldName nf
          bls = CS.lines bs

-- | Read a row of real valued data
readNutRealRow :: Int64 -> Int -> BL.ByteString -> V.Vector Double
readNutRealRow n o = V.fromList . runGet (skip (o * bytesPerReal) >> decoder [])
    where
        o' = fromIntegral (o * bytesPerReal)
        n' = n * bytesPerReal'
        decoder !acc =  do
            !br <- bytesRead
            !be <- isEmpty
            if (br >= (n' + o')) || be
               then return $! reverse acc
               else getDoublebe >>= decoder . (:acc)

-- | Read a row of complex valued data
readNutComplexRow :: Int64 -> Int -> BL.ByteString -> V.Vector (Complex Double)
readNutComplexRow n o = V.fromList . runGet (skip (o * bytesPerComplex) >> decoder)
    where
        o' = fromIntegral (o * bytesPerComplex)
        n' = n * bytesPerComplex'
        decoder =  do
            br <- bytesRead
            be <- isEmpty
            if (br >= (2 * n' + o')) || be
               then return []
               else do !r  <- getDoublebe
                       !c  <- getDoublebe
                       !ds <- decoder
                       return $! ( (r:+c) : ds )

-- vl = [rc (j * nv) | j <- [ 0 .. (np - 1) ]]
-- map V.length vl

-- | Construct NutPlot type
nutPlot :: String -> NutPlotType -> Int -> Int -> [String] -> BL.ByteString -> NutPlot
nutPlot pn nt np nv vn bs
        = NutPlot { nutPlotName  = pn
                  , nutNumVars   = nv
                  , nutNumPoints = np
                  , nutVariables = vn
                  , nutPlotType  = nt
                  , nutData      = M.fromList $ zip vn dt }
    where
      rr = \i -> readNutRealRow (fromIntegral nv) (fromIntegral i) bs
      rc = \i -> readNutComplexRow (fromIntegral nv) (fromIntegral i) bs
      !dt = if nt == NutRealPlot
               then map NutRealWave    . A.toColumns . A.fromRows
                                       $ [rr (j * nv) | j <- [ 0 .. (np - 1) ]]
               else map NutComplexWave . A.toColumns . A.fromRows
                                       $ [rc (j * nv) | j <- [ 0 .. (np - 1) ]]

-- | Get the first NutMeg Element and the Rest
popNutElement :: NutField -> BS.ByteString -> (Maybe BS.ByteString, BS.ByteString)
popNutElement field dat = (elem'', rest)
    where nextField = succ field
          (elem', rest) = BS.breakSubstring (nutFieldName nextField) dat
          elem''        = CS.strip <$> BS.stripPrefix (nutFieldName field) elem'

-- | Extract NutMeg header information
parseNutHeader :: BS.ByteString -> [String]
parseNutHeader = ph NutTitle
    where ph nf bs | BS.null bs = []
                   | (nf == NutTitle) || (nf == NutDate) =
                       let (e, r) = popNutElement nf bs
                        in (CS.unpack . fromJust $ e) : ph (succ nf) r
                   | otherwise = undefined

-- | Extract NutPlot header information.
parseNutPlotHeader :: BS.ByteString -> [String]
parseNutPlotHeader = ph NutPlotname
    where ph nf bs | nf `elem` [ NutPlotname .. NutNoPoints ] =
                       let (e, r) = popNutElement nf bs
                        in (CS.unpack . fromJust $ e) : ph (succ nf) r
                   | otherwise = []

-- | Returns a NutPlot for a given ByteString segment.
parseNutPlot :: BS.ByteString -> NutPlot
parseNutPlot plt = nutPlot pn fl np nv vn dt
    where
      (hdr, dat) = BS.breakSubstring (nutFieldName NutBinData) plt
      hdr'       = parseNutPlotHeader hdr
      pn         = head hdr'
      fl         = nutFlag (hdr' !! 1)
      nv         = read (hdr' !! 2) :: Int
      np         = read (hdr' !! 3) :: Int
      vr         = snd . BS.breakSubstring "\t" . snd
                 . BS.breakSubstring (nutFieldName NutVariables) $ hdr
      vn         = map (CS.unpack . (!!1) . CS.words) $ CS.lines vr
      dt         = BL.fromStrict . fromJust . BS.stripPrefix "Binary:\n" $ dat

-- | Split a nut String at BS.null s
splitNutString :: BS.ByteString -> [NutPlot] -> [NutPlot]
splitNutString bs bl | BS.null bs = bl
                     | otherwise = splitNutString sp (parseNutPlot sl : bl)
                            where Just ps  = BS.stripPrefix pnf bs
                                  (fp, sp) = BS.breakSubstring pnf ps
                                  sl       = BS.concat [pnf, fp]
                                  pnf      = nutFieldName NutPlotname

-- | Returns a NutMeg for given file contents.
parseNutMeg :: BS.ByteString -> NutMeg
parseNutMeg nut = NutMeg { nutTitle = nt
                         , nutDate  = nd
                         , nutPlots = mps }
    where
      (hdr, bdy) = BS.breakSubstring (nutFieldName NutPlotname) nut
      nt         = readNutElement NutTitle hdr
      nd         = readNutElement NutDate hdr
      nps        = splitNutString bdy []
      mps        = M.fromList $ zip (map nutPlotName nps) nps

-- encodeNutPlot :: NutPlot -> ByteString

max_stars_repo_name: AugustUnderground/hnutmeg

max_issues_repo_name: AugustUnderground/hnutmeg

max_forks_repo_name: AugustUnderground/hnutmeg


text:
module Main where

import Criterion.Main (defaultMain, bgroup, bench, nf, nfIO, whnf)
import qualified Numeric.LinearAlgebra as LA
import qualified Numeric.Morpheus as M
import qualified Naive as N

main :: IO ()
main = do
  !ar <- LA.rand 1000 1000
  let !ac = LA.tr ar

  defaultMain [
    bgroup "Activation" [
        bench "sigmoid"            $ nf M.sigmoid ar
        , bench "sigmoid gradient" $ nf M.sigmoidGradient ar
        ]
    , bgroup "Matrix Reduce" [
        bench "sum by columns - row major"          $ nf M.columnSum ar
        , bench "sum by columns - column major"       $ nf M.columnSum ac
        , bench "sum by rows - row major"               $ nf M.rowSum ar
        , bench "sum by rows - column major"          $ nf M.rowSum ac
        , bench "max index in columns - row major"    $ nf M.columnMaxIndex ar
        , bench "max index in columns - column major" $ nf M.columnMaxIndex ac
        , bench "max index in rows - row major"       $ nf M.rowMaxIndex ar
        , bench "max index in rows - column major"    $ nf M.rowMaxIndex ac
        ]
    , bgroup "Matrix Reduce (naive)" [
        bench "max index in columns - row major"    $ nf (N.reduceByColumnsV (fromIntegral . LA.maxIndex) ) ar
        , bench "max index in columns - column major" $ nf (N.reduceByColumnsV (fromIntegral . LA.maxIndex) ) ac
        , bench "max index in rows - row major"         $ nf (N.reduceByRowsV (fromIntegral . LA.maxIndex) ) ar
        , bench "max index in rows - column major"    $ nf (N.reduceByRowsV (fromIntegral . LA.maxIndex) ) ac
        ]
    ]

max_stars_repo_name: Alexander-Ignatyev/morpheus

max_issues_repo_name: aligusnet/morpheus

max_forks_repo_name: aligusnet/morpheus


text:
{-# LANGUAGE NoImplicitPrelude #-}
module Algebra.Additive (
    -- * Class
    C,
    zero,
    (+), (-),
    negate, subtract,

    -- * Complex functions
    sum, sum1,
    sumNestedAssociative,
    sumNestedCommutative,

    -- * Instance definition helpers
    elementAdd, elementSub, elementNeg,
    (<*>.+), (<*>.-), (<*>.-$),

    -- * Instances for atomic types
    propAssociative,
    propCommutative,
    propIdentity,
    propInverse,
  ) where

import qualified Algebra.Laws as Laws

import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  )
import Data.Word (Word, Word8, Word16, Word32, Word64, )

import qualified NumericPrelude.Elementwise as Elem
import Control.Applicative (Applicative(pure, (<*>)), )
import Data.Tuple.HT (fst3, snd3, thd3, )
import qualified Data.List.Match as Match

import qualified Data.Complex as Complex98
import qualified Data.Ratio as Ratio98
import qualified Prelude as P
import Prelude (Integer, Float, Double, fromInteger, )
import NumericPrelude.Base

infixl 6  +, -

{- |
Additive a encapsulates the notion of a commutative group, specified
by the following laws:

@
          a + b === b + a
    (a + b) + c === a + (b + c)
       zero + a === a
   a + negate a === 0
@

Typical examples include integers, dollars, and vectors.

Minimal definition: '+', 'zero', and ('negate' or '(-)')
-}

class C a where
    {-# MINIMAL zero, (+), ((-) | negate) #-}
    -- | zero element of the vector space
    zero     :: a
    -- | add and subtract elements
    (+), (-) :: a -> a -> a
    -- | inverse with respect to '+'
    negate   :: a -> a

    {-# INLINE negate #-}
    negate a = zero - a
    {-# INLINE (-) #-}
    a - b    = a + negate b

{- |
'subtract' is @(-)@ with swapped operand order.
This is the operand order which will be needed in most cases
of partial application.
-}
subtract :: C a => a -> a -> a
subtract = flip (-)

{- |
Sum up all elements of a list.
An empty list yields zero.

This function is inappropriate for number types like Peano.
Maybe we should make 'sum' a method of Additive.
This would also make 'lengthLeft' and 'lengthRight' superfluous.
-}
sum :: (C a) => [a] -> a
sum = foldl (+) zero

{- |
Sum up all elements of a non-empty list.
This avoids including a zero which is useful for types
where no universal zero is available.
-}
sum1 :: (C a) => [a] -> a
sum1 = foldl1 (+)

{- |
Sum the operands in an order,
such that the dependencies are minimized.
Does this have a measurably effect on speed?

Requires associativity.
-}
sumNestedAssociative :: (C a) => [a] -> a
sumNestedAssociative [] = zero
sumNestedAssociative [x] = x
sumNestedAssociative xs = sumNestedAssociative (sum2 xs)

{-
Make sure that the last entries in the list
are equally often part of an addition.
Maybe this can reduce rounding errors.
The list that sum2 computes is a breadth-first-flattened binary tree.

Requires associativity and commutativity.
-}
sumNestedCommutative :: (C a) => [a] -> a
sumNestedCommutative [] = zero
sumNestedCommutative xs@(_:rs) =
   let ys = xs ++ Match.take rs (sum2 ys)
   in  last ys

_sumNestedCommutative :: (C a) => [a] -> a
_sumNestedCommutative [] = zero
_sumNestedCommutative xs@(_:rs) =
   let ys = xs ++ take (length rs) (sum2 ys)
   in  last ys

{-
[a,b,c, a+b,c+(a+b)]
[a,b,c,d, a+b,c+d,(a+b)+(c+d)]
[a,b,c,d,e, a+b,c+d,e+(a+b),(c+d)+e+(a+b)]
[a,b,c,d,e,f, a+b,c+d,e+f,(a+b)+(c+d),(e+f)+((a+b)+(c+d))]
-}

sum2 :: (C a) => [a] -> [a]
sum2 (x:y:rest) = (x+y) : sum2 rest
sum2 xs = xs

{- |
Instead of baking the add operation into the element function,
we could use higher rank types
and pass a generic @uncurry (+)@ to the run function.
We do not do so in order to stay Haskell 98
at least for parts of NumericPrelude.
-}
{-# INLINE elementAdd #-}
elementAdd ::
   (C x) =>
   (v -> x) -> Elem.T (v,v) x
elementAdd f =
   Elem.element (\(x,y) -> f x + f y)

{-# INLINE elementSub #-}
elementSub ::
   (C x) =>
   (v -> x) -> Elem.T (v,v) x
elementSub f =
   Elem.element (\(x,y) -> f x - f y)

{-# INLINE elementNeg #-}
elementNeg ::
   (C x) =>
   (v -> x) -> Elem.T v x
elementNeg f =
   Elem.element (negate . f)

-- like <*>
infixl 4 <*>.+, <*>.-, <*>.-$

{- |
> addPair :: (Additive.C a, Additive.C b) => (a,b) -> (a,b) -> (a,b)
> addPair = Elem.run2 $ Elem.with (,) <*>.+  fst <*>.+  snd
-}
{-# INLINE (<*>.+) #-}
(<*>.+) ::
   (C x) =>
   Elem.T (v,v) (x -> a) -> (v -> x) -> Elem.T (v,v) a
(<*>.+) f acc =
   f <*> elementAdd acc

{-# INLINE (<*>.-) #-}
(<*>.-) ::
   (C x) =>
   Elem.T (v,v) (x -> a) -> (v -> x) -> Elem.T (v,v) a
(<*>.-) f acc =
   f <*> elementSub acc

{-# INLINE (<*>.-$) #-}
(<*>.-$) ::
   (C x) =>
   Elem.T v (x -> a) -> (v -> x) -> Elem.T v a
(<*>.-$) f acc =
   f <*> elementNeg acc

-- * Instances for atomic types

instance C Integer where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Float   where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Double  where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Int     where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Int8    where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Int16   where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Int32   where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Int64   where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Word    where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Word8   where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Word16  where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Word32  where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

instance C Word64  where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = P.fromInteger 0
   negate = P.negate
   (+)    = (P.+)
   (-)    = (P.-)

-- * Instances for composed types

instance (C v0, C v1) => C (v0, v1) where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = (,) zero zero
   (+)    = Elem.run2 $ pure (,) <*>.+  fst <*>.+  snd
   (-)    = Elem.run2 $ pure (,) <*>.-  fst <*>.-  snd
   negate = Elem.run  $ pure (,) <*>.-$ fst <*>.-$ snd

instance (C v0, C v1, C v2) => C (v0, v1, v2) where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero   = (,,) zero zero zero
   (+)    = Elem.run2 $ pure (,,) <*>.+  fst3 <*>.+  snd3 <*>.+  thd3
   (-)    = Elem.run2 $ pure (,,) <*>.-  fst3 <*>.-  snd3 <*>.-  thd3
   negate = Elem.run  $ pure (,,) <*>.-$ fst3 <*>.-$ snd3 <*>.-$ thd3

instance (C v) => C [v] where
   zero   = []
   negate = map negate
   (+) (x:xs) (y:ys) = (+) x y : (+) xs ys
   (+) xs     []     = xs
   (+) []     ys     = ys
   (-) (x:xs) (y:ys) = (-) x y : (-) xs ys
   (-) xs     []     = xs
   (-) []     ys     = negate ys

instance (C v) => C (b -> v) where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero       _ = zero
   (+)    f g x = (+) (f x) (g x)
   (-)    f g x = (-) (f x) (g x)
   negate f   x = negate (f x)

-- * Properties

propAssociative :: (Eq a, C a) => a -> a -> a -> Bool
propCommutative :: (Eq a, C a) => a -> a -> Bool
propIdentity    :: (Eq a, C a) => a -> Bool
propInverse     :: (Eq a, C a) => a -> Bool

propCommutative  =  Laws.commutative (+)
propAssociative  =  Laws.associative (+)
propIdentity     =  Laws.identity (+) zero
propInverse      =  Laws.inverse (+) negate zero

-- legacy

instance (P.Integral a) => C (Ratio98.Ratio a) where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero                =  P.fromInteger 0
   (+)                 =  (P.+)
   (-)                 =  (P.-)
   negate              =  P.negate

instance (P.RealFloat a) => C (Complex98.Complex a) where
   {-# INLINE zero #-}
   {-# INLINE negate #-}
   {-# INLINE (+) #-}
   {-# INLINE (-) #-}
   zero                =  P.fromInteger 0
   (+)                 =  (P.+)
   (-)                 =  (P.-)
   negate              =  P.negate