text:
{-# LANGUAGE BangPatterns, DeriveDataTypeable, RecordWildCards #-}
-- |
-- Module      : Criterion.Analysis
-- Copyright   : (c) 2009-2014 Bryan O'Sullivan
--
-- License     : BSD-style
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : GHC
--
-- Analysis code for benchmarks.

module Criterion.Analysis
    (
      Outliers(..)
    , OutlierEffect(..)
    , OutlierVariance(..)
    , SampleAnalysis(..)
    , analyseSample
    , scale
    , analyseMean
    , countOutliers
    , classifyOutliers
    , noteOutliers
    , outlierVariance
    , resolveAccessors
    , validateAccessors
    , regress
    ) where

-- Temporary: to support pre-AMP GHC 7.8.4:
import Data.Monoid

import Control.Arrow (second)
import Control.Monad (unless, when)
import Control.Monad.Reader (ask)
import Control.Monad.Trans
import Control.Monad.Trans.Except
import Criterion.IO.Printf (note, prolix)
import Criterion.Measurement (secs, threshold)
import Criterion.Monad (Criterion, getGen, getOverhead)
import Criterion.Types
import Data.Int (Int64)
import Data.Maybe (fromJust)
import Statistics.Function (sort)
import Statistics.Quantile (weightedAvg)
import Statistics.Regression (bootstrapRegress, olsRegress)
import Statistics.Resampling (resample)
import Statistics.Sample (mean)
import Statistics.Sample.KernelDensity (kde)
import Statistics.Types (Estimator(..), Sample)
import System.Random.MWC (GenIO)
import qualified Data.List as List
import qualified Data.Map as Map
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import qualified Statistics.Resampling.Bootstrap as B
import Prelude

-- | Classify outliers in a data set, using the boxplot technique.
classifyOutliers :: Sample -> Outliers
classifyOutliers sa = U.foldl' ((. outlier) . mappend) mempty ssa
    where outlier e = Outliers {
                        samplesSeen = 1
                      , lowSevere = if e <= loS && e < hiM then 1 else 0
                      , lowMild = if e > loS && e <= loM then 1 else 0
                      , highMild = if e >= hiM && e < hiS then 1 else 0
                      , highSevere = if e >= hiS && e > loM then 1 else 0
                      }
          !loS = q1 - (iqr * 3)
          !loM = q1 - (iqr * 1.5)
          !hiM = q3 + (iqr * 1.5)
          !hiS = q3 + (iqr * 3)
          q1   = weightedAvg 1 4 ssa
          q3   = weightedAvg 3 4 ssa
          ssa  = sort sa
          iqr  = q3 - q1

-- | Compute the extent to which outliers in the sample data affect
-- the sample mean and standard deviation.
outlierVariance :: B.Estimate  -- ^ Bootstrap estimate of sample mean.
                -> B.Estimate  -- ^ Bootstrap estimate of sample
                               --   standard deviation.
                -> Double      -- ^ Number of original iterations.
                -> OutlierVariance
outlierVariance µ σ a = OutlierVariance effect desc varOutMin
  where
    ( effect, desc ) | varOutMin < 0.01 = (Unaffected, "no")
                     | varOutMin < 0.1  = (Slight,     "slight")
                     | varOutMin < 0.5  = (Moderate,   "moderate")
                     | otherwise        = (Severe,     "severe")
    varOutMin = (minBy varOut 1 (minBy cMax 0 µgMin)) / σb2
    varOut c  = (ac / a) * (σb2 - ac * σg2) where ac = a - c
    σb        = B.estPoint σ
    µa        = B.estPoint µ / a
    µgMin     = µa / 2
    σg        = min (µgMin / 4) (σb / sqrt a)
    σg2       = σg * σg
    σb2       = σb * σb
    minBy f q r = min (f q) (f r)
    cMax x    = fromIntegral (floor (-2 * k0 / (k1 + sqrt det)) :: Int)
      where
        k1    = σb2 - a * σg2 + ad
        k0    = -a * ad
        ad    = a * d
        d     = k * k where k = µa - x
        det   = k1 * k1 - 4 * σg2 * k0

-- | Count the total number of outliers in a sample.
countOutliers :: Outliers -> Int64
countOutliers (Outliers _ a b c d) = a + b + c + d
{-# INLINE countOutliers #-}

-- | Display the mean of a 'Sample', and characterise the outliers
-- present in the sample.
analyseMean :: Sample
            -> Int              -- ^ Number of iterations used to
                                -- compute the sample.
            -> Criterion Double
analyseMean a iters = do
  let µ = mean a
  _ <- note "mean is %s (%d iterations)\n" (secs µ) iters
  noteOutliers . classifyOutliers $ a
  return µ

-- | Multiply the 'Estimate's in an analysis by the given value, using
-- 'B.scale'.
scale :: Double                 -- ^ Value to multiply by.
      -> SampleAnalysis -> SampleAnalysis
scale f s@SampleAnalysis{..} = s {
                                 anMean = B.scale f anMean
                               , anStdDev = B.scale f anStdDev
                               }

-- | Perform an analysis of a measurement.
analyseSample :: Int            -- ^ Experiment number.
              -> String         -- ^ Experiment name.
              -> V.Vector Measured -- ^ Sample data.
              -> ExceptT String Criterion Report
analyseSample i name meas = do
  Config{..} <- ask
  overhead <- lift getOverhead
  let ests      = [Mean,StdDev]
      -- The use of filter here throws away very-low-quality
      -- measurements when bootstrapping the mean and standard
      -- deviations.  Without this, the numbers look nonsensical when
      -- very brief actions are measured.
      stime     = measure (measTime . rescale) .
                  G.filter ((>= threshold) . measTime) . G.map fixTime .
                  G.tail $ meas
      fixTime m = m { measTime = measTime m - overhead / 2 }
      n         = G.length meas
      s         = G.length stime
  _ <- lift $ prolix "bootstrapping with %d of %d samples (%d%%)\n"
              s n ((s * 100) `quot` n)
  gen <- lift getGen
  rs <- mapM (\(ps,r) -> regress gen ps r meas) $
        ((["iters"],"time"):regressions)
  resamps <- liftIO $ resample gen ests resamples stime
  let [estMean,estStdDev] = B.bootstrapBCA confInterval stime ests resamps
      ov = outlierVariance estMean estStdDev (fromIntegral n)
      an = SampleAnalysis {
               anRegress    = rs
             , anOverhead   = overhead
             , anMean       = estMean
             , anStdDev     = estStdDev
             , anOutlierVar = ov
             }
  return Report {
      reportNumber   = i
    , reportName     = name
    , reportKeys     = measureKeys
    , reportMeasured = meas
    , reportAnalysis = an
    , reportOutliers = classifyOutliers stime
    , reportKDEs     = [uncurry (KDE "time") (kde 128 stime)]
    }

-- | Regress the given predictors against the responder.
--
-- Errors may be returned under various circumstances, such as invalid
-- names or lack of needed data.
--
-- See 'olsRegress' for details of the regression performed.
regress :: GenIO
        -> [String]             -- ^ Predictor names.
        -> String               -- ^ Responder name.
        -> V.Vector Measured
        -> ExceptT String Criterion Regression
regress gen predNames respName meas = do
  when (G.null meas) $
    throwE "no measurements"
  accs <- ExceptT . return $ validateAccessors predNames respName
  let unmeasured = [n | (n, Nothing) <- map (second ($ G.head meas)) accs]
  unless (null unmeasured) $
    throwE $ "no data available for " ++ renderNames unmeasured
  let (r:ps)      = map ((`measure` meas) . (fromJust .) . snd) accs
  Config{..} <- ask
  (coeffs,r2) <- liftIO $
                 bootstrapRegress gen resamples confInterval olsRegress ps r
  return Regression {
      regResponder = respName
    , regCoeffs    = Map.fromList (zip (predNames ++ ["y"]) (G.toList coeffs))
    , regRSquare   = r2
    }

singleton :: [a] -> Bool
singleton [_] = True
singleton _   = False

-- | Given a list of accessor names (see 'measureKeys'), return either
-- a mapping from accessor name to function or an error message if
-- any names are wrong.
resolveAccessors :: [String]
                 -> Either String [(String, Measured -> Maybe Double)]
resolveAccessors names =
  case unresolved of
    [] -> Right [(n, a) | (n, Just (a,_)) <- accessors]
    _  -> Left $ "unknown metric " ++ renderNames unresolved
  where
    unresolved = [n | (n, Nothing) <- accessors]
    accessors = flip map names $ \n -> (n, Map.lookup n measureAccessors)

-- | Given predictor and responder names, do some basic validation,
-- then hand back the relevant accessors.
validateAccessors :: [String]   -- ^ Predictor names.
                  -> String     -- ^ Responder name.
                  -> Either String [(String, Measured -> Maybe Double)]
validateAccessors predNames respName = do
  when (null predNames) $
    Left "no predictors specified"
  let names = respName:predNames
      dups = map head . filter (not . singleton) .
             List.group . List.sort $ names
  unless (null dups) $
    Left $ "duplicated metric " ++ renderNames dups
  resolveAccessors names

renderNames :: [String] -> String
renderNames = List.intercalate ", " . map show

-- | Display a report of the 'Outliers' present in a 'Sample'.
noteOutliers :: Outliers -> Criterion ()
noteOutliers o = do
  let frac n = (100::Double) * fromIntegral n / fromIntegral (samplesSeen o)
      check :: Int64 -> Double -> String -> Criterion ()
      check k t d = when (frac k > t) $
                    note "  %d (%.1g%%) %s\n" k (frac k) d
      outCount = countOutliers o
  when (outCount > 0) $ do
    _ <- note "found %d outliers among %d samples (%.1g%%)\n"
         outCount (samplesSeen o) (frac outCount)
    check (lowSevere o) 0 "low severe"
    check (lowMild o) 1 "low mild"
    check (highMild o) 1 "high mild"
    check (highSevere o) 0 "high severe"

max_stars_repo_name: Shimuuar/criterion

max_issues_repo_name: Shimuuar/criterion

max_forks_repo_name: Shimuuar/criterion


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

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

import ArrayFire

spec :: Spec
spec =
  describe "Array tests" $ do
    it "Should perform Array tests" $ do
      (1 + 1) `shouldBe` 2
    it "Should fail to create 0 dimension arrays" $ do
      let arr = mkArray @Int [0,0,0,0] [1..]
      evaluate arr `shouldThrow` anyException
    it "Should fail to create 0 length arrays" $ do
      let arr = mkArray @Int [0,0,0,1] []
      evaluate arr `shouldThrow` anyException
    it "Should fail to create 0 length arrays w/ 0 dimensions" $ do
      let arr = mkArray @Int [0,0,0,0] []
      evaluate arr `shouldThrow` anyException
    it "Should create a column vector" $ do
      let arr = mkArray @Int [9,1,1,1] (repeat 9)
      isColumn arr `shouldBe` True
    it "Should create a row vector" $ do
      let arr = mkArray @Int [1,9,1,1] (repeat 9)
      isRow arr `shouldBe` True
    it "Should create a vector" $ do
      let arr = mkArray @Int [9,1,1,1] (repeat 9)
      isVector arr `shouldBe` True
    it "Should create a vector" $ do
      let arr = mkArray @Int [1,9,1,1] (repeat 9)
      isVector arr `shouldBe` True
    it "Should copy an array" $ do
      let arr = mkArray @Int [9,9,1,1] (repeat 9)
      let newArray = copyArray arr
      newArray `shouldBe` arr
    it "Should modify manual eval flag" $ do
      setManualEvalFlag False
      (`shouldBe` False) =<< getManualEvalFlag
    it "Should return the number of elements" $ do
      let arr = mkArray @Int [9,9,1,1] [1..]
      getElements arr `shouldBe` 81
--    it "Should give an empty array" $ do
--      let arr = mkArray @Int [-1,1,1,1] []
--      getElements arr `shouldBe` 0
--      isEmpty arr `shouldBe` True
    it "Should create a scalar array" $ do
      let arr = mkArray @Int [1] [1]
      isScalar arr `shouldBe` True
    it "Should get number of dims specified" $ do
      let arr = mkArray @Int [1,1,1,1] [1]
      getNumDims arr `shouldBe` 1
      let arr = mkArray @Int [2,3,4,5] [1..]
      getNumDims arr `shouldBe` 4
      let arr = mkArray @Int [2,3,4] [1..]
      getNumDims arr `shouldBe` 3
    it "Should get value of dims specified" $ do
      let arr = mkArray @Int [2,3,4,5] (repeat 1)
      getDims arr `shouldBe` (2,3,4,5)

    it "Should test Sparsity" $ do
      let arr = mkArray @Double [2,2,1,1] (repeat 1)
      isSparse arr `shouldBe` False

    it "Should make a Bit array" $ do
      let arr = mkArray @CBool [2,2] [1,1,1,1]
      isBool arr `shouldBe` True

    it "Should make an integer array" $ do
      let arr = mkArray @Int [2,2] (repeat 1)
      isInteger arr `shouldBe` True

    it "Should make a Floating array" $ do
      let arr = mkArray @Double [2,2] (repeat 1)
      isFloating arr `shouldBe` True
      let arr = mkArray @CBool [2,2] (repeat 1)
      isFloating arr `shouldBe` False

    it "Should make a Complex array" $ do
      let arr = mkArray @(Complex Double) [2,2] (repeat 1)
      isComplex arr `shouldBe` True
      isReal arr `shouldBe` False

    it "Should make a Real array" $ do
      let arr = mkArray @Double [2,2] (repeat 1)
      isReal arr `shouldBe` True
      isComplex arr `shouldBe` False

    it "Should make a Double precision array" $ do
      let arr = mkArray @Double [2,2] (repeat 1)
      isDouble arr `shouldBe` True
      isSingle arr `shouldBe` False

    it "Should make a Single precision array" $ do
      let arr = mkArray @Float [2,2] (repeat 1)
      isDouble arr `shouldBe` False
      isSingle arr `shouldBe` True

    it "Should make a Real floating array" $ do
      let arr = mkArray @Float [2,2] (repeat 1)
      isRealFloating arr `shouldBe` True
      let arr = mkArray @Double [2,2] (repeat 1)
      isRealFloating arr `shouldBe` True

    it "Should get reference count" $ do
      let arr1 = mkArray @Float [2,2] (repeat 1)
          arr2 = retainArray arr1
          arr3 = retainArray arr2
      getDataRefCount arr3 `shouldBe` 3

    it "Should convert an array to a list" $ do
      let arr = mkArray @Double [30,30] (repeat 1)
      toList arr `shouldBe` Prelude.replicate (30 * 30) 1

      let arr = mkArray @Float [10,10] (repeat (5.5))
      toList arr `shouldBe` Prelude.replicate 100 5.5

      let arr = mkArray @CBool [4] [1,1,0,1]
      toList arr `shouldBe` [1,1,0,1]

      let arr = mkArray @Int16 [10] [1..]
      toList arr `shouldBe` [1..10]

      let arr = mkArray @Int32 [100] [1..100]
      toList arr `shouldBe` [1..100]

      let arr = mkArray @Int64 [100] [1..100]
      toList arr `shouldBe` [1..100]

      let arr = mkArray @Int [100] [1..100]
      toList arr `shouldBe` [1..100]

      let arr = mkArray @(Complex Float) [1] [1 :+ 1]
      toList arr `shouldBe` [1 :+ 1]

      let arr = mkArray @(Complex Double) [1] [1 :+ 1]
      toList arr `shouldBe` [1 :+ 1]

      let arr = mkArray @Word16 [10] [1..10]
      toList arr `shouldBe` [1..10]

      let arr = mkArray @Word32 [10] [1..10]
      toList arr `shouldBe` [1..10]

      let arr = mkArray @Word64 [10] [1..10]
      toList arr `shouldBe` [1..10]

      let arr = mkArray @Word [10] [1..10]
      toList arr `shouldBe` [1..10]

ext: hs

max_stars_repo_name: danielkroeni/arrayfire-haskell

max_issues_repo_name: danielkroeni/arrayfire-haskell

max_forks_repo_name: danielkroeni/arrayfire-haskell


text:
-- | This module provides generic functionality to deal with ensembles in
-- statistical mechanics.

module StatisticalMechanics.Ensemble where

import Numeric.Log

import Statistics.Probability

-- | The state probability functions provide conversion from some types @a@
-- into non-normalized probabilities. For "real" applications, using the
-- @logProbability@ function is preferred. This functions allows for easy
-- abstraction when types @a@ are given as fractions of some actual value (say:
-- deka-cal), or are discretized.
--
-- The returned values are not normalized, because we do not now the total
-- evidence @Z@ until integration over all states has happened -- which is not
-- feasible in a number of problems.
--
-- TODO replace @()@ with temperature and results with non-normalized @P@ or
-- @LogP@, depending. At some point we want to have type-level physical
-- quantities, hence the need for the second type.

class StateProbability a where
  -- | Given a temperature and a state "energy", return the corresponding
  -- non-normalized probability.
  stateProbability
    ∷ Double
    -- ^ this is @k*T@
    → a
    -- ^ the energy (or discretized energy)
    → Probability NotNormalized Double
    -- ^ probability of being in state @a@, but only proportional up to @1/Z@.
  stateLogProbability
    ∷ Double
    -- ^ this is @1/(k * T)@
    → a
    -- ^ the energy (or discretized energy)
    → Log (Probability NotNormalized Double)
    -- ^ resulting probability

instance StateProbability Double where
  stateProbability kT x = Prob . exp . negate $ x/kT
  {-# Inline stateProbability #-}
  --stateLogProbability kT x = Exp . log . Prob . exp . negate $ x/kT
  stateLogProbability kT x = Exp . Prob . negate $ x/kT
  {-# Inline stateLogProbability #-}

max_stars_repo_name: choener/SciBaseTypes

max_issues_repo_name: choener/SciBaseTypes

max_forks_repo_name: choener/SciBaseTypes


text:
import Criterion.Main
import Numeric.LinearAlgebra

n = 10 :: Int
main = defaultMain [
         bench "matrix creation" $ w (\_ -> fromLists sampleList :: Matrix Double),
         bench "random matrix generation" $ whnfIO (randn n n),
         bench "matrix multiplication" $ w (\_ -> sampleDiag <> sampleDiag),
         bench "matrix slicing" $ whnf (\n -> sampleDiag ?? (Take (n - 1), Take (n - 1))) n,

         bench "scalar multiplication" $ whnf (\n -> sampleDiag * (fromIntegral n)) n,

         bench "matrix-vector addition" $ w (\_ -> sampleDiag + sampleDiag),

         bench "transpose" $ w (\_ -> tr sampleDiag)
       ]
  where
    w f = f `whnf` 0
    sampleList = replicate n $ replicate n 1
    sampleDiag = diagl [1..fromIntegral n]

max_stars_repo_name: DataHaskell/numeric-libs-benchmarks

max_issues_repo_name: mdibaiee/numeric-libs-overview

max_forks_repo_name: DataHaskell/numeric-libs-benchmarks


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ViewPatterns     #-}
{-|
Module      : Grenade.Layers.Internal.Hmatrix
Description : Faster implementations of HMatrix functions
Maintainer  : Theo Charalambous
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Internal.Hmatrix (
    reshapeMatrix
  ) where

import           Numeric.LinearAlgebra
import           Numeric.LinearAlgebra.Devel

-- | Reshape a row major matrix
reshapeMatrix :: (Element t, Num t, Container Vector t) => Int -> Int -> Matrix t -> Matrix t
reshapeMatrix r c m@(size->(r', c'))
    | r * c == r' * c' = matrixFromVector RowMajor r c $ flatten m
    | otherwise        = error $ "can't reshape matrix of shape dim = " ++ show (r', c') ++ " to matrix of shape " ++ show (r, c)
{-# INLINE reshapeMatrix #-}

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE GADTs, RecordWildCards, ScopedTypeVariables #-}
module Data.Eigen.SparseMatrix.Mutable (
    -- * Mutable SparseMatrix
    IOSparseMatrix(..),
    IOSparseMatrixXf,
    IOSparseMatrixXd,
    IOSparseMatrixXcf,
    IOSparseMatrixXcd,
    new,
    reserve,
    -- * Matrix properties
    rows,
    cols,
    innerSize,
    outerSize,
    nonZeros,
    -- * Matrix compression
    compressed,
    compress,
    uncompress,
    -- * Accessing matrix data
    read,
    write,
    setZero,
    setIdentity,
    -- * Changing matrix shape
    resize,
    conservativeResize
) where

import Prelude hiding (read)
import Data.Complex
import Foreign.C.String
import Foreign.C.Types
import Foreign.ForeignPtr
import Foreign.Marshal.Alloc
import Foreign.Ptr
import Foreign.Storable
import qualified Foreign.Concurrent as FC
import qualified Data.Eigen.Internal as I

-- | Mutable version of sparse matrix. See `Data.Eigen.SparseMatrix.SparseMatrix` for details about matrix layout.
data IOSparseMatrix a b where
    IOSparseMatrix :: I.Elem a b => !(ForeignPtr (I.CSparseMatrix a b)) -> IOSparseMatrix a b

-- | Alias for single precision mutable matrix
type IOSparseMatrixXf = IOSparseMatrix Float CFloat
-- | Alias for double precision mutable matrix
type IOSparseMatrixXd = IOSparseMatrix Double CDouble
-- | Alias for single previsiom mutable matrix of complex numbers
type IOSparseMatrixXcf = IOSparseMatrix (Complex Float) (I.CComplex CFloat)
-- | Alias for double prevision mutable matrix of complex numbers
type IOSparseMatrixXcd = IOSparseMatrix (Complex Double) (I.CComplex CDouble)

-- | Creates new matrix with the given size @rows x cols@
new :: I.Elem a b => Int -> Int -> IO (IOSparseMatrix a b)
new rows cols = alloca $ \pm -> do
    I.call $ I.sparse_new (I.cast rows) (I.cast cols) pm
    m <- peek pm
    fm <- FC.newForeignPtr m $ I.call $ I.sparse_free m
    return $! IOSparseMatrix fm

-- | Returns the number of rows of the matrix
rows :: I.Elem a b => IOSparseMatrix a b -> IO Int
rows = _prop I.sparse_rows (return . I.cast)

-- | Returns the number of columns of the matrix
cols :: I.Elem a b => IOSparseMatrix a b  -> IO Int
cols = _prop I.sparse_cols (return . I.cast)

-- | Returns the number of rows (resp. columns) of the matrix if the storage order column major (resp. row major)
innerSize :: I.Elem a b => IOSparseMatrix a b  -> IO Int
innerSize = _prop I.sparse_innerSize (return . I.cast)

-- | Returns the number of columns (resp. rows) of the matrix if the storage order column major (resp. row major)
outerSize :: I.Elem a b => IOSparseMatrix a b  -> IO Int
outerSize = _prop I.sparse_outerSize (return . I.cast)

-- | Returns whether this matrix is in compressed form.
compressed :: I.Elem a b => IOSparseMatrix a b -> IO Bool
compressed = _prop I.sparse_isCompressed (return . (==1))

-- | Turns the matrix into the compressed format.
compress :: I.Elem a b => IOSparseMatrix a b -> IO ()
compress = _inplace I.sparse_compressInplace

-- | Turns the matrix into the uncompressed mode.
uncompress :: I.Elem a b => IOSparseMatrix a b -> IO ()
uncompress = _inplace I.sparse_uncompressInplace

-- | Reads the value of the matrix at position @i@, @j@.
-- This function returns @Scalar(0)@ if the element is an explicit zero.
read :: I.Elem a b => IOSparseMatrix a b -> Int -> Int -> IO a
read (IOSparseMatrix fm) row col = withForeignPtr fm $ \m -> alloca $ \px -> do
    I.call $ I.sparse_coeff m (I.cast row) (I.cast col) px
    I.cast <$> peek px

{- | Writes the value of the matrix at position @i@, @j@.
    This function turns the matrix into a non compressed form if that was not the case.

    This is a @O(log(nnz_j))@ operation (binary search) plus the cost of element insertion if the element does not already exist.

    Cost of element insertion is sorted insertion in O(1) if the elements of each inner vector are inserted in increasing inner index order, and in @O(nnz_j)@ for a random insertion.
-}
write :: I.Elem a b => IOSparseMatrix a b -> Int -> Int -> a -> IO ()
write (IOSparseMatrix fm) row col x = withForeignPtr fm $ \m -> alloca $ \px -> do
    I.call $ I.sparse_coeffRef m (I.cast row) (I.cast col) px
    peek px >>= (`poke` I.cast x)

-- | Sets the matrix to the identity matrix
setIdentity :: I.Elem a b => IOSparseMatrix a b -> IO ()
setIdentity = _inplace I.sparse_setIdentity

-- | Removes all non zeros but keep allocated memory
setZero :: I.Elem a b => IOSparseMatrix a b -> IO ()
setZero = _inplace I.sparse_setZero

-- | The number of non zero coefficients
nonZeros :: I.Elem a b => IOSparseMatrix a b -> IO Int
nonZeros = _prop I.sparse_nonZeros (return . I.cast)

-- | Preallocates space for non zeros. The matrix must be in compressed mode.
reserve :: I.Elem a b => IOSparseMatrix a b -> Int -> IO ()
reserve m s = _inplace (\p -> I.sparse_reserve p (I.cast s)) m

-- | Resizes the matrix to a rows x cols matrix and initializes it to zero.
resize :: I.Elem a b => IOSparseMatrix a b -> Int -> Int -> IO ()
resize m rows cols = _inplace (\p -> I.sparse_resize p (I.cast rows) (I.cast cols)) m

-- | Resizes the matrix to a rows x cols matrix leaving old values untouched.
conservativeResize :: I.Elem a b => IOSparseMatrix a b -> Int -> Int -> IO ()
conservativeResize m rows cols = _inplace (\p -> I.sparse_conservativeResize p (I.cast rows) (I.cast cols)) m

_inplace :: I.Elem a b => (Ptr (I.CSparseMatrix a b) -> IO CString) -> IOSparseMatrix a b -> IO ()
_inplace f (IOSparseMatrix fm) = withForeignPtr fm $ \m -> I.call $ f m

_prop :: Storable c => (I.CSparseMatrixPtr a b -> Ptr c -> IO CString) -> (c -> IO d) -> IOSparseMatrix a b -> IO d
_prop f g (IOSparseMatrix fp) =
    withForeignPtr fp $ \p ->
        alloca $ \pq -> do
            I.call (f p pq)
            peek pq >>= g

max_stars_repo_name: osidorkin/haskell-eigen

max_issues_repo_name: osidorkin/haskell-eigen

max_forks_repo_name: osidorkin/haskell-eigen


text:
{-# LANGUAGE OverloadedStrings #-}

import           Chapter04.LinearAlgebra.ListBasedVectorProperties (linearAlgebraTests)
import           StatisticsTests                                   (statsSpecs)
import           Test.Tasty

main :: IO ()
main = do
  statsSpecs' <- statsSpecs
  let tests = testGroup "all tests" [ statsSpecs', linearAlgebraTests ]
  defaultMain tests

currentTests :: IO ()
currentTests = statsSpecs >>= \x ->
  defaultMain $ testGroup "current development modules" [ x ]


text:
{-# LANGUAGE NamedFieldPuns, TemplateHaskell #-}

module School.App.Test.FileTransformer
( fileTransformerTest ) where

import Conduit ((.|), await, liftIO, runConduit, sinkList)
import qualified Data.ByteString.Lazy as BL
import Control.Applicative (liftA2)
import Data.Default.Class (def)
import Data.Function (on)
import Numeric.LinearAlgebra ((?), Matrix, R)
import School.App.FileTransformer
import School.FileIO.AppIO (runAppIO)
import School.FileIO.FileApp (fileApp)
import School.FileIO.FileType (FileType(..), toExtension)
import School.FileIO.FileHeader (FileHeader(..))
import School.FileIO.Source (source)
import School.TestUtils (assertRight, dummyMatrix)
import School.Types.DataType (DataType(..))
import School.Types.Error (Error)
import School.Utils.Either (isLeft, isRight)
import System.Directory (removeFile)
import Test.Tasty (TestTree)
import Test.Tasty.QuickCheck
import Test.Tasty.TH
import Test.QuickCheck.Monadic (PropertyM, assert, monadicIO, run)

inFileName :: FileType -> FilePath
inFileName fType = "test/data/matrix3x3."
                ++ toExtension fType

inFileNameTwice :: DataType -> FileType -> FilePath
inFileNameTwice DBL64B fType = "test/data/matrix3x3Twice."
                            ++ toExtension fType
inFileNameTwice INT32B fType = "test/data/matrix3x3Twice_INT32B."
                            ++ toExtension fType
inFileNameTwice INT08B fType = "test/data/matrix3x3Twice_INT08B."
                            ++ toExtension fType

header3x3 :: FileHeader
header3x3 = FileHeader { dataType = DBL64B
                         , cols = 3
                         , rows = 3
                         }

testFile :: FileType -> String
testFile fType = "test." ++ toExtension fType

fileEq :: FilePath -> FilePath -> IO Bool
fileEq = liftA2 (==) `on` BL.readFile

readMat :: FileHeader
        -> FilePath
        -> PropertyM IO (Either Error
                                (Maybe (Matrix R)))
readMat header path = run . runAppIO . runConduit $
  source header path .| await

readDoubleList :: FileHeader
               -> FilePath
               -> PropertyM IO (Either Error
                                       [Double])
readDoubleList header path = run . runAppIO . runConduit $
  source header path .| sinkList

prop_copy :: FileType -> Property
prop_copy fType = monadicIO $ do
  let options = def { inFileOpt = inFileName fType
                    , inFileTypeOpt = Just fType
                    , outFileOpt = testFile fType
                    , outFileTypeOpt = Just fType
                    }
  result <- liftIO . runAppIO $ fileApp options
  assert $ isRight result
  equal <- liftIO $ fileEq (inFileName fType) (testFile fType)
  liftIO $ removeFile (testFile fType)
  assert equal

prop_copy_guess_filetypes :: FileType -> Property
prop_copy_guess_filetypes fType = monadicIO $ do
  let options = def { inFileOpt = (inFileName fType)
                    , outFileOpt = (testFile fType)
                    }
  result <-liftIO . runAppIO $ fileApp options
  assert $ isRight result
  equal <- liftIO $ fileEq (inFileName fType) (testFile fType)
  liftIO $ removeFile (testFile fType)
  assert equal

prop_copy_add_extension :: FileType -> Property
prop_copy_add_extension fType = monadicIO $ do
  let options = def { inFileOpt = inFileName fType
                    , inFileTypeOpt = Just fType
                    , outFileOpt = "test"
                    , outFileTypeOpt = Just fType
                    }
  result <- liftIO . runAppIO $ fileApp options
  assert $ isRight result
  equal <- liftIO $ fileEq (inFileName fType) (testFile fType)
  liftIO $ removeFile (testFile fType)
  assert equal

prop_skip_rows :: FileType -> Property
prop_skip_rows fType = monadicIO $ do
  let outHeader = header3x3 { rows = 1 }
  let options = def { inFileOpt = inFileName fType
                    , inFileTypeOpt = Just fType
                    , outFileOpt = testFile fType
                    , outFileTypeOpt = Just fType
                    , skipRowsOpt = Just 2
                    }
  result <- liftIO . runAppIO $ fileApp options
  assert $ isRight result
  let check = Just $ dummyMatrix 3 3 ? [2]
  readRes <- readMat outHeader (testFile fType)
  liftIO $ removeFile (testFile fType)
  assertRight (== check) readRes

prop_skip_too_many_rows :: FileType -> Property
prop_skip_too_many_rows fType = monadicIO $ do
  let options = def { inFileOpt = inFileName fType
                    , inFileTypeOpt = Just fType
                    , skipRowsOpt = Just 4
                    }
  result <- liftIO . runAppIO $ fileApp options
  assert $ isLeft result

prop_limit_extent :: FileType -> Property
prop_limit_extent fType = monadicIO $ do
  let options = def { inFileOpt = inFileName fType
                    , inFileTypeOpt = Just fType
                    , nRowsOpt = Just 2
                    , outFileOpt = testFile fType
                    }
  result <-liftIO . runAppIO $ fileApp options
  assert $ isRight result
  let check = Just $ dummyMatrix 3 3 ? [0, 1]
  let outHeader = header3x3 { rows = 2 }
  readRes <- readMat outHeader (testFile fType)
  liftIO $ removeFile (testFile fType)
  assertRight (== check) readRes

prop_skip_and_limit :: FileType -> Property
prop_skip_and_limit fType = monadicIO $ do
  let options = def { inFileOpt = inFileName fType
                    , inFileTypeOpt = Just fType
                    , nRowsOpt = Just 1
                    , skipRowsOpt = Just 1
                    , outFileOpt = testFile fType
                    }
  result <-liftIO . runAppIO $ fileApp options
  assert $ isRight result
  let check = Just $ dummyMatrix 3 3 ? [1]
  let outHeader = header3x3 { rows = 1 }
  readRes <- readMat outHeader (testFile fType)
  liftIO $ removeFile (testFile fType)
  assertRight (== check) readRes

prop_limit_too_many_rows :: FileType -> Property
prop_limit_too_many_rows fType = monadicIO $ do
  let options = def { inFileOpt = inFileName fType
                    , inFileTypeOpt = Just fType
                    , nRowsOpt = Just 4
                    }
  result <- liftIO . runAppIO $ fileApp options
  assert $ isLeft result

prop_limit_and_skip_too_many :: FileType -> Property
prop_limit_and_skip_too_many fType = monadicIO $ do
  let options = def { inFileOpt = inFileName fType
                    , inFileTypeOpt = Just SM
                    , nRowsOpt = Just 3
                    , skipRowsOpt = Just 1
                    }
  result <- liftIO . runAppIO $ fileApp options
  assert $ isLeft result

inFileType :: FileType -> DataType -> FilePath
inFileType fType DBL64B = "test/data/matrix3x3."
                       ++ toExtension fType
inFileType fType INT32B = "test/data/matrix3x3_INT32B."
                       ++ toExtension fType
inFileType fType INT08B = "test/data/matrix3x3_INT08B."
                       ++ toExtension fType

legalConversion :: DataType -> DataType -> Bool
legalConversion INT08B _ = True
legalConversion INT32B INT32B = True
legalConversion INT32B DBL64B = True
legalConversion DBL64B DBL64B = True
legalConversion _ _ = False

prop_convert_data_type :: FileType -> DataType -> DataType -> Property
prop_convert_data_type fType tIn tOut = monadicIO $ do
  let options = def { inFileOpt = inFileType fType tIn
                    , inDataTypeOpt = tIn
                    , inFileTypeOpt = Just fType
                    , outDataTypeOpt = Just tOut
                    , outFileOpt = testFile fType
                    , outFileTypeOpt = Just fType
                    }
  result <- liftIO . runAppIO $ fileApp options
  if legalConversion tIn tOut
    then do
      assert $ isRight result
      let inHeader = header3x3 { dataType = tIn }
      let outHeader = header3x3 { dataType = tOut }
      inRes <- readDoubleList inHeader (inFileType fType tIn)
      outRes <- readDoubleList outHeader (testFile fType)
      liftIO $ removeFile (testFile fType)
      assert $ inRes == outRes
    else assert $ isLeft result

prop_convert_file_type :: FileType -> FileType -> DataType -> Property
prop_convert_file_type fTypeIn fTypeOut dType = monadicIO $ do
  let options = def { inFileOpt = inFileType fTypeIn dType
                    , inDataTypeOpt = dType
                    , inFileTypeOpt = Just fTypeIn
                    , outFileOpt = testFile fTypeOut
                    , outFileTypeOpt = Just fTypeOut
                    }
  result <- liftIO . runAppIO $ fileApp options
  assert $ isRight result
  let header = header3x3 { dataType = dType }
  inRes <- readDoubleList header (inFileType fTypeIn dType)
  outRes <- readDoubleList header (testFile fTypeOut)
  liftIO $ removeFile (testFile fTypeOut)
  assert $ inRes == outRes

prop_convert_both_types :: FileType -> FileType -> DataType -> DataType -> Property
prop_convert_both_types fTypeIn fTypeOut tIn tOut = monadicIO $ do
  let options = def { inFileOpt = inFileType fTypeIn tIn
                    , inDataTypeOpt = tIn
                    , inFileTypeOpt = Just fTypeIn
                    , outDataTypeOpt = Just tOut
                    , outFileOpt = testFile fTypeOut
                    , outFileTypeOpt = Just fTypeOut
                    }
  result <- liftIO . runAppIO $ fileApp options
  if legalConversion tIn tOut
    then do
      assert $ isRight result
      let inHeader = header3x3 { dataType = tIn }
      let outHeader = header3x3 { dataType = tOut }
      inRes <- readDoubleList inHeader (inFileType fTypeIn tIn)
      outRes <- readDoubleList outHeader (testFile fTypeOut)
      liftIO $ removeFile (testFile fTypeOut)
      assert $ inRes == outRes
    else assert $ isLeft result

prop_copy_twice :: FileType -> Property
prop_copy_twice fType = monadicIO $ do
  let options = def { inFileOpt = inFileNameTwice DBL64B fType
                    , inFileTypeOpt = Just fType
                    , outFileOpt = testFile fType
                    , outFileTypeOpt = Just fType
                    }
  result <- liftIO . runAppIO $ fileApp options
  assert $ isRight result
  let header = header3x3 { rows = 6 }
  inRes <- readDoubleList header (inFileNameTwice DBL64B fType)
  outRes <- readDoubleList header (testFile fType)
  liftIO $ removeFile (testFile fType)
  assert $ inRes == outRes

prop_convert_data_twice :: FileType -> DataType -> DataType -> Property
prop_convert_data_twice fType tIn tOut = monadicIO $ do
  let options = def { inFileOpt = inFileNameTwice tIn fType
                    , inDataTypeOpt = tIn
                    , inFileTypeOpt = Just fType
                    , outDataTypeOpt = Just tOut
                    , outFileOpt = testFile fType
                    , outFileTypeOpt = Just fType
                    }
  result <- liftIO . runAppIO $ fileApp options
  if legalConversion tIn tOut
    then do
      assert $ isRight result
      let inHeader = header3x3 { dataType = tIn }
      let outHeader = header3x3 { dataType = tOut }
      inRes <- readDoubleList inHeader (inFileNameTwice tIn fType)
      outRes <- readDoubleList outHeader (testFile fType)
      liftIO $ removeFile (testFile fType)
      assert $ inRes == outRes
    else do
      liftIO $ removeFile (testFile fType)
      assert $ isLeft result

fileTransformerTest :: TestTree
fileTransformerTest = $(testGroupGenerator)

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
{-# LANGUAGE TypeFamilies, GADTs, EmptyDataDecls, RankNTypes, EmptyCase #-}
{-# LANGUAGE DataKinds, ConstraintKinds, PolyKinds #-}
{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts #-}
{-# LANGUAGE TypeOperators, UndecidableInstances, ScopedTypeVariables #-}
{-# LANGUAGE FunctionalDependencies, DeriveFunctor, StandaloneDeriving #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# OPTIONS_GHC -Wno-redundant-constraints #-}

module Language.Grappa.GrappaInternals where

import Data.List
import Data.Typeable
import Data.Functor.Const
import Data.Proxy (Proxy(..))
import GHC.Exts (Constraint, IsList(..))

import Data.Aeson
import Data.Aeson.Types (Parser)
import Data.Vector (Vector)
import qualified Data.Vector as V
import qualified Numeric.LinearAlgebra as M

import Language.Grappa.Distribution

--
-- * Representing Grappa Types
--

-- | A type is "atomic" iff it is not an 'ADT'
type family IsAtomic a :: Bool where
  IsAtomic (ADT g) = 'False
  IsAtomic a = 'True

-- | The Grappa type of distributions. How it is represented depends on the
-- interpretation being used, so we do not instantiate it to an actual Haskell
-- type.
data Dist a deriving Typeable

-- | Is a type considered to be a base type in Grappa?
type family IsBaseType a where
  IsBaseType (a -> b) = 'False
  IsBaseType (ADT a) = 'False
  IsBaseType (Dist a) = 'False
  IsBaseType a = 'True

-- | An object-level representation of lists of Grappa types
data GrappaTypeListRepr (as :: [*]) where
  GrappaTypeListNil :: GrappaTypeListRepr '[]
  GrappaTypeListCons :: GrappaTypeRepr a -> GrappaTypeListRepr as ->
                        GrappaTypeListRepr (a ': as)

-- | An object-level representation of an application of a type constructor to 0
-- or more Grappa types
data GrappaTypeAppRepr (a :: k) where
  GrappaTypeAppBase :: Typeable a => GrappaTypeAppRepr a
  GrappaTypeAppApply :: GrappaTypeAppRepr (f :: * -> k) ->
                        GrappaTypeRepr a ->
                        GrappaTypeAppRepr (f a)

-- | An object-level representation of Grappa types
data GrappaTypeRepr (a :: *) where
  GrappaBaseType :: (IsAtomic a ~ 'True, IsBaseType a ~ 'True) =>
                    GrappaTypeAppRepr a -> GrappaTypeRepr a
  -- ^ A "base type" imported directly from Haskell, applied to 0 or more Grappa
  -- types; it must not be equal to any other Grappa type construct
  GrappaADTType :: GrappaADT adt => GrappaTypeAppRepr adt ->
                   GrappaTypeRepr (ADT adt)
  -- ^ A Grappa ADT type
  GrappaTupleType ::
    GrappaADT (TupleF as) =>
    GrappaTypeListRepr as -> GrappaTypeRepr (ADT (TupleF as))
  -- ^ A Grappa tuple type; note that we require 'GrappaADT', since the
  -- instance for 'TupleF' requires a class instance of 'GrappaTypeList' for
  -- @as@, not just a 'GrappaTypeListRepr'
  GrappaDistType :: GrappaTypeRepr a -> GrappaTypeRepr (Dist a)
  -- ^ A Grappa distribution type
  GrappaArrowType :: GrappaTypeRepr a -> GrappaTypeRepr b ->
                     GrappaTypeRepr (a -> b)
  -- ^ A Grappa function type

instance Show (GrappaTypeRepr a) where
  show _ = "(FIXME HERE: write a Show instance for GrappaTypeRepr!)"

-- | A typeclass indicating that a type is a valid Grappa type
class GrappaType a where
  grappaTypeRepr :: GrappaTypeRepr a

-- | A typeclass indicating that a list of types are all valid Grappa types. We
-- represent it as a type family as well as a type constraint, below, as the
-- type family helps GHC pull out all the 'GrappaType' instances for each elem
type family GrappaTypeListFam (as :: [*]) :: Constraint where
  GrappaTypeListFam '[] = ()
  -- To help GHC not have to unroll as many type applications, we put some
  -- helper cases to skip forward multiple types in a type list
  GrappaTypeListFam (a ': '[]) = GrappaType a
  GrappaTypeListFam (a ': b ': '[]) = (GrappaType a, GrappaType b)
  GrappaTypeListFam (a ': b ': c ': '[]) =
    (GrappaType a, GrappaType b, GrappaType c)
  GrappaTypeListFam (a ': b ': c ': d ': '[]) =
    (GrappaType a, GrappaType b, GrappaType c, GrappaType d)
  GrappaTypeListFam (a ': b ': c ': d ': e ': rest) =
    (GrappaType a, GrappaType b, GrappaType c, GrappaType d,
     GrappaType e, GrappaTypeList rest)

class (IsTypeList as, GrappaTypeListFam as) => GrappaTypeList as where
  grappaTypeListRepr :: GrappaTypeListRepr as

-- Instances for lists of types
instance GrappaTypeList '[] where
  grappaTypeListRepr = GrappaTypeListNil
instance GrappaType a => GrappaTypeList '[a] where
  grappaTypeListRepr = GrappaTypeListCons grappaTypeRepr GrappaTypeListNil
instance (GrappaType a, GrappaType b) => GrappaTypeList '[a,b] where
  grappaTypeListRepr =
    GrappaTypeListCons grappaTypeRepr
    (GrappaTypeListCons grappaTypeRepr GrappaTypeListNil)
instance (GrappaType a, GrappaType b, GrappaType c) =>
         GrappaTypeList '[a,b,c] where
  grappaTypeListRepr =
    GrappaTypeListCons grappaTypeRepr
    (GrappaTypeListCons grappaTypeRepr
     (GrappaTypeListCons grappaTypeRepr GrappaTypeListNil))
instance (GrappaType a, GrappaType b, GrappaType c, GrappaType d) =>
         GrappaTypeList '[a,b,c,d] where
  grappaTypeListRepr =
    GrappaTypeListCons grappaTypeRepr
    (GrappaTypeListCons grappaTypeRepr
     (GrappaTypeListCons grappaTypeRepr
      (GrappaTypeListCons grappaTypeRepr GrappaTypeListNil)))
instance (GrappaType a, GrappaType b, GrappaType c,
          GrappaType d, GrappaType e, GrappaTypeList rest) =>
         GrappaTypeList (a ': b ': c ': d ': e ': rest) where
  grappaTypeListRepr =
    GrappaTypeListCons grappaTypeRepr
    (GrappaTypeListCons grappaTypeRepr
     (GrappaTypeListCons grappaTypeRepr
      (GrappaTypeListCons grappaTypeRepr
       (GrappaTypeListCons grappaTypeRepr grappaTypeListRepr))))

-- Instances for representing the core Grappa type constructs
instance GrappaTypeList as => GrappaType (ADT (TupleF as)) where
  grappaTypeRepr = GrappaTupleType grappaTypeListRepr
instance GrappaType a => GrappaType (Dist a) where
  grappaTypeRepr = GrappaDistType grappaTypeRepr
instance (GrappaType a, GrappaType b) => GrappaType (a -> b) where
  grappaTypeRepr = GrappaArrowType grappaTypeRepr grappaTypeRepr
instance GrappaType a => GrappaType (Vector a) where
  grappaTypeRepr =
    GrappaBaseType $ GrappaTypeAppApply GrappaTypeAppBase grappaTypeRepr

-- Instances for representing base types; remember to add more instances here
-- when we want to support more Haskell types in Grappa!
instance GrappaType Bool where
  grappaTypeRepr = GrappaBaseType GrappaTypeAppBase
instance GrappaType Int where
  grappaTypeRepr = GrappaBaseType GrappaTypeAppBase
instance GrappaType R where
  grappaTypeRepr = GrappaBaseType GrappaTypeAppBase
instance GrappaType Prob where
  grappaTypeRepr = GrappaBaseType GrappaTypeAppBase

-- Instances for real and probability vectors and matrices
instance GrappaType RVector where
  grappaTypeRepr = GrappaBaseType GrappaTypeAppBase
instance GrappaType RMatrix where
  grappaTypeRepr = GrappaBaseType GrappaTypeAppBase
instance GrappaType ProbVector where
  grappaTypeRepr = GrappaBaseType GrappaTypeAppBase
instance GrappaType ProbMatrix where
  grappaTypeRepr = GrappaBaseType GrappaTypeAppBase

-- | Test if a 'GrappaTypeRepr' is equal to 'R'
matchGrappaRType :: GrappaTypeRepr a -> Maybe (a :~: R)
matchGrappaRType (GrappaBaseType GrappaTypeAppBase) = eqT
matchGrappaRType _ = Nothing

-- | Test if a 'GrappaTypeRepr' is equal to 'Int'
matchGrappaIntType :: GrappaTypeRepr a -> Maybe (a :~: Int)
matchGrappaIntType (GrappaBaseType GrappaTypeAppBase) = eqT
matchGrappaIntType _ = Nothing

--
-- * Grappa ADTs
--

-- | Helper type: the identity
newtype Id a = Id {unId :: a} deriving (Functor, Show, Num, Floating, Fractional)

instance Applicative Id where
  pure = Id
  f <*> x = Id $ unId f $ unId x

-- FIXME: This approach to modeling ADTs does not let the type arguments change
-- in recursive uses of the ADT

-- | Builds an ADT type from an ADT type name
newtype ADT (adt :: (* -> *) -> * -> *) =
  ADT { unADT :: adt Id (ADT adt) }

-- | Class for mapping functions over ADTs
class TraversableADT (adt :: (* -> *) -> * -> *) where
  traverseADT :: Applicative m => (forall a. f a -> m (g a)) ->
                 adt f r -> m (adt g r)

-- | Map a function over an 'ADT'
mapADT :: TraversableADT adt => (forall a. f a -> g a) -> adt f r -> adt g r
mapADT f = unId . traverseADT (Id . f)

-- | Fold a function over all the constructor arguments of an ADT. Note that
-- this is /not/ the same as, say, a list fold, since here the "rest" argument
-- in a list is treated as a single argument, i.e., we do not recurse.
foldrADT :: TraversableADT adt =>
            (forall a. f a -> r) -> (r -> r -> r) -> r -> adt f (ADT adt) -> r
foldrADT resF f x = foldr f x . ctorArgsADT resF

-- | A function that can test if a given ADT uses a specific constructor
newtype CtorMatcher adt =
  CtorMatcher { applyCtorMatcher :: forall f. adt f (ADT adt) -> Bool }

-- | Typeclass of ADTs where we can enumerate all the constructors, by building
-- a list of ADT elements with 'Proxy's as arguments
class ReifyCtorsADT adt where
  reifyCtorsADT :: [adt Proxy (ADT adt)]

-- | Convert all the arguments of the constructor of an ADT into a list
ctorArgsADT :: TraversableADT adt => (forall a. f a -> r) ->
               adt f (ADT adt) -> [r]
ctorArgsADT resF = getConst . traverseADT (Const . (: []) . resF)

-- | Class for "valid" Grappa ADT types
class (TraversableADT adt, ReifyCtorsADT adt) =>
      GrappaADT (adt :: (* -> *) -> * -> *) where
  -- | Stronger form of 'TraversableADT', where we know that the argument types
  -- are themselves valid Grappa types
  gtraverseADT :: Applicative m =>
                  (forall a. GrappaType a => f a -> m (g a)) ->
                  adt f (ADT adt) -> m (adt g (ADT adt))

-- | Map a function on well-typed objects over an 'ADT'
gmapADT :: GrappaADT adt => (forall a. GrappaType a => f a -> g a) ->
           adt f (ADT adt) -> adt g (ADT adt)
gmapADT f = unId . gtraverseADT (Id . f)

-- | Helper type class for 'Show'ing 'ADT's
class ShowADT (adt :: (* -> *) -> * -> *) where
  showADT :: ADT adt -> String

instance ShowADT adt => Show (ADT adt) where
  show = showADT

-- NOTE: we don't want the generalized Show instance for an ADT, but usually
-- want adt-specific printing...
--
-- deriving instance (Show (adt f (ADT adt f))) => Show (ADT adt f)

--
-- * Heterogeneous Lists
--

-- | An @'HList' f [t0, t1, ..., tn]@ is a heterogeneous list of elements of
-- type @f t0@, @f t1@, etc.
data HList (f :: * -> *) (ts  :: [*]) where
  HListNil  :: HList f '[]
  HListCons :: f a -> HList f as -> HList f (a ': as)

-- | Convert a Grappa tuple to an 'HList'
tupleToHList :: TupleF ts f r -> HList f ts
tupleToHList Tuple0 = HListNil
tupleToHList tup@(Tuple1 _)           = consTupleToHList tup
tupleToHList tup@(Tuple2 _ _)         = consTupleToHList tup
tupleToHList tup@(Tuple3 _ _ _)       = consTupleToHList tup
tupleToHList tup@(Tuple4 _ _ _ _)     = consTupleToHList tup
tupleToHList tup@(TupleN _ _ _ _ _ _) = consTupleToHList tup

-- | Helper for 'tupleToHList'
consTupleToHList :: TupleF (t ': ts) f r -> HList f (t ': ts)
consTupleToHList tup =
  HListCons (tupleHead tup) (tupleToHList (tupleTail tup))

-- | Convert an 'HList' to a Grappa tuple
hListToTuple :: HList f ts -> TupleF ts f r
hListToTuple HListNil = Tuple0
hListToTuple (HListCons x xs) = tupleCons x (hListToTuple xs)

-- | Convert an 'HList' to a monomorphic list
hListToList :: (forall a. f a -> b) -> HList f ts -> [b]
hListToList _ HListNil = []
hListToList g (HListCons x xs) = g x : hListToList g xs

-- | Map a function over an 'HList'
mapHList :: (forall a. f a -> g a) -> HList f as -> HList g as
mapHList _ HListNil = HListNil
mapHList f (HListCons a as) = HListCons (f a) (mapHList f as)

-- | Traverse an 'HList'
traverseHList :: Applicative m => (forall a. f a -> m (g a)) ->
                 HList f ts -> m (HList g ts)
traverseHList _ HListNil = pure HListNil
traverseHList f (HListCons x xs) =
  pure HListCons <*> f x <*> traverseHList f xs

-- | Traverse an 'HList' in a monomorphic way
traverseHList' :: Applicative m => (forall a. f a -> m b) -> HList f ts -> m [b]
traverseHList' _ HListNil         = pure []
traverseHList' f (HListCons x xs) =
  pure (:) <*> f x <*> traverseHList' f xs

-- | Fold an 'HList' in a monomorphic way
foldHList :: (forall a. f a -> b -> b) -> HList f ts -> b -> b
foldHList _ HListNil acc         = acc
foldHList g (HListCons x xs) acc = g x $ foldHList g xs acc

-- | Compute the length of an 'HList'
hListLength :: HList f ts -> Int
hListLength l = foldHList (const (+ 1)) l 0

-- | Build an 'HList' of 'TypeListElem' proofs that each type in @as@ is an
-- element of @as@
buildHListElems :: HList f as -> HList (TypeListElem as) as
buildHListElems HListNil = HListNil
buildHListElems (HListCons _ rest) =
  HListCons TypeListElem_Base (mapHList TypeListElem_Cons
                               (buildHListElems rest))

{-
type family Concat (l1 :: [*]) (l2 :: [*]) where
  Concat '[] l = l
  Concat (x:xs) l = x : Concat xs l

hListConcat :: HList f as -> HList f bs -> HList f (Concat as bs)
hListConcat HListNil l = l
hListConcat (HListCons x xs) l = HListCons x (hListConcat xs l)
-}

--
-- * Tuple ADTs
--

-- | A generic tuple type as an ADT
data TupleF (ts :: [*]) f (r :: *) where
  Tuple0 :: TupleF '[] f r
  -- | Note: one-tuples are needed for recursive 'TupleN' constructions
  Tuple1 :: f a -> TupleF '[a] f r
  Tuple2 :: f a -> f b -> TupleF '[a, b] f r
  Tuple3 :: f a -> f b -> f c -> TupleF '[a, b, c] f r
  Tuple4 :: f a -> f b -> f c -> f d -> TupleF '[a, b, c, d] f r
  TupleN :: f a -> f b -> f c -> f d -> f e -> TupleF rest f r ->
            TupleF (a ': b ': c ': d ': e ': rest) f r

type family TupleConstraintF (cf :: * -> Constraint) (f :: * -> *)
     (ts :: [*]) :: Constraint where
  TupleConstraintF cf f '[] = ()
  TupleConstraintF cf f '[a] = (cf (f a))
  TupleConstraintF cf f '[a,b] = (cf (f a), cf (f b))
  TupleConstraintF cf f '[a,b,c] = (cf (f a), cf (f b), cf (f c))
  TupleConstraintF cf f '[a,b,c,d] = (cf (f a), cf (f b), cf (f c), cf (f d))
  TupleConstraintF cf f (a ': b ': c ': d ': e ': ts) =
    (cf (f a), cf (f b), cf (f c), cf (f d), cf (f e),
     TupleConstraintF cf f ts)

deriving instance TupleConstraintF Eq f ts => Eq (TupleF ts f r)

instance (IsTypeList ts, TupleConstraintF FromJSON f ts) =>
         FromJSON (TupleF ts f r) where
  parseJSON = withArray "Grappa tuple" (helper typeListProxy) where
    helper :: TupleConstraintF FromJSON f ts' => TupleF ts' Proxy r ->
              Vector Value -> Parser (TupleF ts' f r)
    helper Tuple0 _ = return Tuple0
    helper (Tuple1 _) vals = Tuple1 <$> parseJSON (vals V.! 0)
    helper (Tuple2 _ _) vals =
      Tuple2 <$> parseJSON (vals V.! 0) <*> parseJSON (vals V.! 1)
    helper (Tuple3 _ _ _) vals =
      Tuple3 <$> parseJSON (vals V.! 0) <*> parseJSON (vals V.! 1)
      <*> parseJSON (vals V.! 2)
    helper (Tuple4 _ _ _ _) vals =
      Tuple4 <$> parseJSON (vals V.! 0) <*> parseJSON (vals V.! 1)
      <*> parseJSON (vals V.! 2) <*> parseJSON (vals V.! 3)
    helper (TupleN _ _ _ _ _ tup) vals =
      TupleN <$> parseJSON (vals V.! 0) <*> parseJSON (vals V.! 1)
      <*> parseJSON (vals V.! 2) <*> parseJSON (vals V.! 3)
      <*> parseJSON (vals V.! 4) <*> helper tup (V.drop 5 vals)

-- | This says that @ts@ is a well-formed list of types, that we can reflect on
class IsTypeList ts where
  typeListProxy :: TupleF ts Proxy r

instance IsTypeList '[] where
  typeListProxy = Tuple0

instance IsTypeList '[a] where
  typeListProxy = Tuple1 Proxy

instance IsTypeList '[a, b] where
  typeListProxy = Tuple2 Proxy Proxy

instance IsTypeList '[a, b, c] where
  typeListProxy = Tuple3 Proxy Proxy Proxy

instance IsTypeList '[a, b, c, d] where
  typeListProxy = Tuple4 Proxy Proxy Proxy Proxy

instance IsTypeList rest => IsTypeList (a ': b ': c ': d ': e ': rest) where
  typeListProxy = TupleN Proxy Proxy Proxy Proxy Proxy typeListProxy

-- | "Proofs" that @t@ is in the list @ts@ of types
data TypeListElem ts t where
  TypeListElem_Base :: TypeListElem (t ': ts) t
  TypeListElem_Cons :: TypeListElem ts t -> TypeListElem (u ': ts) t

-- | Add an extra element to a tuple
tupleCons :: f t -> TupleF ts f r1 -> TupleF (t ': ts) f r2
tupleCons t Tuple0 = Tuple1 t
tupleCons t (Tuple1 a) = Tuple2 t a
tupleCons t (Tuple2 a b) = Tuple3 t a b
tupleCons t (Tuple3 a b c) = Tuple4 t a b c
tupleCons t (Tuple4 a b c d) = TupleN t a b c d Tuple0
tupleCons t (TupleN a b c d e rest) = TupleN t a b c d $ tupleCons e rest

-- | Get the first element of a tuple
tupleHead :: TupleF (t ': ts) f r -> f t
tupleHead (Tuple1 a) = a
tupleHead (Tuple2 a _) = a
tupleHead (Tuple3 a _ _) = a
tupleHead (Tuple4 a _ _ _) = a
tupleHead (TupleN a _ _ _ _ _) = a

-- | Remove the first element of a tuple
tupleTail :: TupleF (t ': ts) f r1 -> TupleF ts f r2
tupleTail (Tuple1 _) = Tuple0
tupleTail (Tuple2 _ b) = Tuple1 b
tupleTail (Tuple3 _ b c) = Tuple2 b c
tupleTail (Tuple4 _ b c d) = Tuple3 b c d
tupleTail (TupleN _ b c d e rest) =
  tupleCons b $ tupleCons c $ tupleCons d $ tupleCons e rest

-- | Project an element of a tuple
projectTuple :: TypeListElem ts t -> TupleF ts f r -> f t
projectTuple TypeListElem_Base tup = tupleHead tup
projectTuple (TypeListElem_Cons elemPf) tup =
  projectTuple elemPf $ tupleTail tup

-- | Project an element of an 'HList'
projectHList :: HList f as -> TypeListElem as a -> f a
projectHList (HListCons t _) TypeListElem_Base = t
projectHList (HListCons _ ts) (TypeListElem_Cons pf) = projectHList ts pf
projectHList HListNil tp_elem = case tp_elem of { }

-- Need a TraversableADT instance for each ADT type
instance TraversableADT (TupleF as) where
  traverseADT _ Tuple0 = pure Tuple0
  traverseADT f (Tuple1 a) = pure Tuple1 <*> f a
  traverseADT f (Tuple2 a b) = pure Tuple2 <*> f a <*> f b
  traverseADT f (Tuple3 a b c) = pure Tuple3 <*> f a <*> f b <*> f c
  traverseADT f (Tuple4 a b c d) = pure Tuple4 <*> f a <*> f b <*> f c <*> f d
  traverseADT f (TupleN a b c d e rest) =
    pure TupleN <*> f a <*> f b <*> f c <*> f d <*> f e <*> traverseADT f rest

untraverseTuple :: (Functor g, IsTypeList bs) =>
                   (forall a. g (f a) -> h a) ->
                   g (TupleF bs f r) -> TupleF bs h r
untraverseTuple = helper typeListProxy
  where
    helper :: Functor g => TupleF bs proxy r ->
              (forall a. g (f a) -> h a) -> g (TupleF bs f r) -> TupleF bs h r
    helper Tuple0 = \_ _ -> Tuple0
    helper tup@(Tuple1 _) = helperStep tup
    helper tup@(Tuple2 _ _) = helperStep tup
    helper tup@(Tuple3 _ _ _) = helperStep tup
    helper tup@(Tuple4 _ _ _ _) = helperStep tup
    helper tup@(TupleN _ _ _ _ _ _) = helperStep tup

    helperStep :: Functor g => TupleF (b ': bs) proxy r ->
                  (forall a. g (f a) -> h a) ->
                  g (TupleF (b ': bs) f r) -> TupleF (b ': bs) h r
    helperStep ts f gtup =
      tupleCons (f $ fmap tupleHead gtup)
      (helper (tupleTail ts) f $ fmap tupleTail gtup)

-- | Build a tuple from a polymorphic function for each element
buildTuple :: IsTypeList bs => (forall a. f a) -> TupleF bs f r
buildTuple f = untraverseTuple (\_ -> f) (Id typeListProxy)

-- Need a ReifyCtorsADT instance for each ADT type
instance IsTypeList ts => ReifyCtorsADT (TupleF ts) where
  reifyCtorsADT = [buildTuple Proxy]

-- | Map a binary function over tuples
mapTuple2 :: (forall a. f a -> g a -> h a) ->
             TupleF ts f r -> TupleF ts g r -> TupleF ts h r
mapTuple2 _ Tuple0 Tuple0 = Tuple0
mapTuple2 f (Tuple1 a1) (Tuple1 a2) = Tuple1 (f a1 a2)
mapTuple2 f (Tuple2 a1 b1) (Tuple2 a2 b2) = Tuple2 (f a1 a2) (f b1 b2)
mapTuple2 f (Tuple3 a1 b1 c1) (Tuple3 a2 b2 c2) =
  Tuple3 (f a1 a2) (f b1 b2) (f c1 c2)
mapTuple2 f (Tuple4 a1 b1 c1 d1) (Tuple4 a2 b2 c2 d2) =
  Tuple4 (f a1 a2) (f b1 b2) (f c1 c2) (f d1 d2)
mapTuple2 f (TupleN a1 b1 c1 d1 e1 rest1) (TupleN a2 b2 c2 d2 e2 rest2) =
  TupleN (f a1 a2) (f b1 b2) (f c1 c2) (f d1 d2) (f e1 e2) $
  mapTuple2 f rest1 rest2

-- Need a GrappaADT instance for each ADT type
instance GrappaTypeList as => GrappaADT (TupleF as) where
  gtraverseADT = helper where
    helper :: (GrappaTypeList bs, Applicative m) =>
              (forall a. GrappaType a => f a -> m (g a)) ->
              TupleF bs f r -> m (TupleF bs g r)
    helper _ Tuple0 = pure Tuple0
    helper f (Tuple1 a) = pure Tuple1 <*> f a
    helper f (Tuple2 a b) = pure Tuple2 <*> f a <*> f b
    helper f (Tuple3 a b c) = pure Tuple3 <*> f a <*> f b <*> f c
    helper f (Tuple4 a b c d) = pure Tuple4 <*> f a <*> f b <*> f c <*> f d
    helper f (TupleN a b c d e rest) =
      pure TupleN <*> f a <*> f b <*> f c <*> f d <*> f e <*> helper f rest

-- | Type synonym for Grappa tuples
type GTuple ts = ADT (TupleF ts)

-- | Map a type function over a list of types
type family MapF (f :: * -> *) (xs :: [*]) where
  MapF f '[] = '[]
  MapF f (x ': xs) = f x ': MapF f xs

-- | Defined type class for mappping a constraint function over a list; note
-- that we special-case small-sized lists, to help GHC do less unfolding
type family MapC (f :: * -> Constraint) (ts :: [*]) :: Constraint where
  MapC f '[] = ()
  MapC f '[a] = f a
  MapC f '[a,b] = (f a, f b)
  MapC f '[a,b,c] = (f a, f b, f c)
  MapC f '[a,b,c,d] = (f a, f b, f c, f d)
  MapC f (a ': b ': c ': d ': e ': rest) =
    (f a, f b, f c, f d, f e, MapC f rest)

instance MapC Show ts => ShowADT (TupleF ts) where
  showADT (ADT tup_body) = helper tup_body where
    helper :: forall ts' r. MapC Show ts' => TupleF ts' Id r -> String
    helper Tuple0 = "()"
    helper (Tuple1 (Id a)) = "(" ++ show a ++ ")"
    helper (Tuple2 (Id a) (Id b)) =
      "(" ++ show a ++ "," ++ show b ++ ")"
    helper (Tuple3 (Id a) (Id b) (Id c)) =
      "(" ++ show a ++ "," ++ show b ++ "," ++ show c ++ ")"
    helper (Tuple4 (Id a) (Id b) (Id c) (Id d)) =
      "(" ++ show a ++ "," ++ show b ++ "," ++ show c ++ "," ++ show d ++ ")"
    helper (TupleN (Id a) (Id b) (Id c) (Id d) (Id e) rest) =
      "(" ++ show a ++ "," ++ show b ++ "," ++ show c ++ "," ++ show d
      ++ "," ++ show e ++ "," ++ helper rest ++ ")"

--
-- * Boolean ADT
--

-- | The list type functor
data BoolF f r = TrueF | FalseF
  deriving (Eq, Show, Typeable)

-- Need a GrappaTypeRepr instance
instance GrappaType (ADT BoolF) where
  grappaTypeRepr = GrappaADTType GrappaTypeAppBase

-- Need a TraversableADT instance for each ADT type
instance TraversableADT BoolF where
  traverseADT _ TrueF  = pure TrueF
  traverseADT _ FalseF = pure FalseF

-- Need a ReifyCtorsADT instance for each ADT type
instance ReifyCtorsADT BoolF where
  reifyCtorsADT = [TrueF, FalseF]

-- Need a GrappaADT instance for each ADT type
instance GrappaADT BoolF where
  gtraverseADT _ TrueF  = pure TrueF
  gtraverseADT _ FalseF = pure FalseF

-- Also need CtorMatchers for each constructor
ctorMatcher__TrueF :: CtorMatcher BoolF
ctorMatcher__TrueF =
  CtorMatcher (\ x ->
               case x of
                 TrueF -> True
                 _ -> False)

ctorMatcher__FalseF :: CtorMatcher BoolF
ctorMatcher__FalseF =
  CtorMatcher (\ x ->
               case x of
                 FalseF -> True
                 _ -> False)

-- | Type synonym for lists, as Grappa ADTs
type GBool = ADT BoolF

fromHaskellBool :: Bool -> BoolF f r
fromHaskellBool True  = TrueF
fromHaskellBool False = FalseF

ifF :: BoolF f a -> t -> t -> t
ifF TrueF  t _ = t
ifF FalseF _ e = e

--
-- * List ADT
--

-- | The list type functor
data ListF a f r = Nil | Cons (f a) (f r)

-- Need a GrappaTypeRepr instance
instance GrappaType a => GrappaType (ADT (ListF a)) where
  grappaTypeRepr =
    GrappaADTType (GrappaTypeAppApply GrappaTypeAppBase grappaTypeRepr)

--deriving instance (Show (f a), Show (f r)) => Show (ListF a f r)

deriving instance (Eq (f a), Eq (f r)) => Eq (ListF a f r)

instance (FromJSON (f a), FromJSON (f r)) => FromJSON (ListF a f r) where
  parseJSON = withArray "Grappa list" $ \vals ->
    if V.null vals then return Nil else
      Cons <$> parseJSON (V.head vals) <*> parseJSON (Array $ V.tail vals)

-- Need a TraversableADT instance for each ADT type
instance TraversableADT (ListF a) where
  traverseADT _ Nil = pure Nil
  traverseADT f (Cons x xs) = pure Cons <*> f x <*> f xs

-- Need a ReifyCtorsADT instance for each ADT type
instance ReifyCtorsADT (ListF a) where
  reifyCtorsADT = [Nil, Cons Proxy Proxy]

-- Need a GrappaADT instance for each ADT type
instance GrappaType a => GrappaADT (ListF a) where
  gtraverseADT _ Nil = pure Nil
  gtraverseADT f (Cons x xs) = pure Cons <*> f x <*> f xs

-- Also need CtorMatchers for each constructor
ctorMatcher__Nil :: CtorMatcher (ListF a)
ctorMatcher__Nil =
  CtorMatcher (\ x ->
               case x of
                 Nil -> True
                 _ -> False)

ctorMatcher__Cons :: CtorMatcher (ListF a)
ctorMatcher__Cons =
  CtorMatcher (\ x ->
               case x of
                 Cons _ _ -> True
                 _ -> False)

-- | Type synonym for lists, as Grappa ADTs
type GList a = ADT (ListF a)

instance Show a => ShowADT (ListF a) where
  showADT = show . toList

listRepeatF :: a -> ADT (ListF a)
listRepeatF f = x
  where x = ADT (Cons (pure f) (pure x))

listIterateF :: a -> (a -> a) -> ADT (ListF a)
listIterateF x f = ADT (Cons (pure x) (pure (listIterateF (f x) f)))

enumFromByF :: (Num a) => a -> a -> ADT (ListF a)
enumFromByF from by = listIterateF from (+by)

enumFromToByF :: (Ord a, Num a) => a -> a -> a -> ADT (ListF a)
enumFromToByF from to by = go from
  where go n
          | n > to    = ADT Nil
          | otherwise = ADT (Cons (pure n) (pure (go (n + by))))

fromHaskellList :: [a] -> GList a
fromHaskellList = fromList

toHaskellList :: GList a -> [a]
toHaskellList = toList

-- | Convert a Grappa list with a non-'Id' functor to a Haskell list
toHaskellListF :: (f (ADT (ListF a)) -> ListF a f (ADT (ListF a))) ->
                  ListF a f (ADT (ListF a)) -> [f a]
toHaskellListF _ Nil = []
toHaskellListF unwrap (Cons x xs) = x : toHaskellListF unwrap (unwrap xs)

-- | Convert a Haskell list with a non-'Id' functor to a Grappa list
fromHaskellListF :: (ListF a f (ADT (ListF a)) -> f (ADT (ListF a))) ->
                    [f a] -> ListF a f (ADT (ListF a))
fromHaskellListF _ [] = Nil
fromHaskellListF wrap (x : xs) = Cons x (wrap $ fromHaskellListF wrap xs)

instance IsList (ADT (ListF a)) where
  type Item (ADT (ListF a)) = a
  fromList [] = ADT Nil
  fromList (x:xs) = ADT (Cons (Id x) (Id (fromList xs)))
  toList (ADT Nil) = []
  toList (ADT (Cons (Id x) (Id xs))) = x : toList xs

--
-- * Grappa-Specific Printing Methods
--

class GrappaShow t where
  grappaShow :: t -> String

instance GrappaShow Double where
  grappaShow = show

instance GrappaShow Prob where
  grappaShow = show

instance GrappaShow Int where
  grappaShow = show

instance GrappaShow a => GrappaShow (Id a) where
  grappaShow (Id x) = grappaShow x

instance GrappaShow (adt Id (ADT adt)) => GrappaShow (ADT adt) where
  grappaShow (ADT x) = grappaShow x

instance GrappaShow a => GrappaShow (Vector a) where
  grappaShow xs =
    "[|" ++ concat (intersperse ", " (map grappaShow $ V.toList xs)) ++ "|]"

instance GrappaShow RVector where
  grappaShow (RVector xs) =
    "[|" ++ concat (intersperse ", " (map grappaShow $ M.toList xs)) ++ "|]"

instance GrappaShow ProbVector where
  grappaShow (ProbVector xs) =
    "[|" ++ concat (intersperse ", "
                    (map (grappaShow . logRToProb) $ M.toList xs)) ++ "|]"

instance GrappaShow RMatrix where
  grappaShow m =
    "[R|" ++ concat (intersperse "\n" $ map grappaShow (toRowsM m)) ++ "|]"

instance GrappaShow ProbMatrix where
  grappaShow m =
    "[P|" ++ concat (intersperse "\n" $ map grappaShow (toRowsPM m)) ++ "|]"

instance MapC GrappaShow (MapF f ts) => GrappaShow (TupleF ts f r) where
  grappaShow Tuple0 = "()"
  grappaShow (Tuple1 a) = "(" ++ grappaShow a ++ ")"
  grappaShow (Tuple2 a b) = "(" ++ grappaShow a ++ "," ++ grappaShow b ++ ")"
  grappaShow (Tuple3 a b c) =
    "(" ++ grappaShow a ++ "," ++ grappaShow b ++ "," ++ grappaShow c ++ ")"
  grappaShow (Tuple4 a b c d) =
    "(" ++ grappaShow a ++ "," ++ grappaShow b ++ "," ++ grappaShow c ++
    "," ++ grappaShow d ++ ")"
  grappaShow (TupleN a b c d e _) =
    "(" ++ grappaShow a ++ "," ++ grappaShow b ++ "," ++ grappaShow c ++
    "," ++ grappaShow d ++ "," ++ grappaShow e ++ ", ...)"

instance (GrappaShow (f t), GrappaShowListContents (f r)) => GrappaShow (ListF t f r) where
  grappaShow Nil = "[]"
  grappaShow lst = "[" ++ go (showListContents lst) ++ "]"
    where go [] = ""
          go [x] = x
          go (x:xs) = x ++ "," ++ go xs

class GrappaShowListContents t where
  showListContents :: t -> [String]

instance (GrappaShow (f t), GrappaShowListContents (f r)) => GrappaShowListContents (ListF t f r) where
  showListContents Nil = []
  showListContents (Cons x xs) =
    grappaShow x : showListContents xs

instance GrappaShowListContents t => GrappaShowListContents (Id t) where
  showListContents (Id x) = showListContents x

max_stars_repo_name: kquick/grappa

max_issues_repo_name: kquick/grappa

max_forks_repo_name: kquick/grappa


text:
{-# LANGUAGE Haskell2010, TemplateHaskell #-}

module Numeric.Matrix.Sugar (
        iMatrix, dMatrix, cMatrix
    ) where

import Numeric.Matrix
import Language.Haskell.TH
import Language.Haskell.TH.Quote

import Data.Data
import Data.Complex
import Data.Ratio

type ReadM a = String -> Matrix a

iMatrix = QuasiQuoter
            (quoter (read :: ReadM Integer))
            undefined undefined undefined

dMatrix = QuasiQuoter
            (quoter (read :: ReadM Double))
            undefined undefined undefined

cMatrix = QuasiQuoter
            (quoter (read :: ReadM (Complex Double)))
            undefined undefined undefined

quoter :: (Data e, MatrixElement e)
       => (String -> Matrix e) -> String -> Q Exp

quoter read str = do
    let qExp = dataToExpQ (const Nothing) (toList (read str))
    [e| fromList $qExp |]

max_stars_repo_name: phadej/bed-and-breakfast

max_issues_repo_name: phadej/bed-and-breakfast

max_forks_repo_name: phadej/bed-and-breakfast


text:
module Statistics.Information.Utils.Matrix where

import Data.Matrix
import qualified Data.Vector as Vector
import Statistics.Information.Utils.List

rows :: Matrix a -> [Vector.Vector a]
rows xs = [getRow i xs | i <- [1..nrows xs]]

columns :: Matrix a -> [Vector.Vector a]
columns xs = [getCol i xs | i <- [1..ncols xs]]

rowExtend :: Matrix a -> Int -> Matrix a
rowExtend xs n =
  if remainder > 0 then
    remRows <-> Prelude.foldr (<->) xs (repeatN (completions - 1) xs)
  else
    Prelude.foldr (<->) xs (repeatN (completions - 1) xs)
  where
    completions = n `quot` nrows xs
    remRows = fromLists [Vector.toList $ getRow i xs |
                         i <- [1..remainder]]
    remainder = n `rem` nrows xs

residualMatrix :: Matrix a -> Int -> Matrix a
xsiLeft xs i = submatrix 1 (nrows xs) 1 i xs
xsiRight xs i = submatrix 1 (nrows xs) (i+1) (ncols xs) xs
residualMatrix xs i | i == (ncols xs) = xsiLeft xs i
residualMatrix xs 1 = xsiRight xs 1
residualMatrix xs i = xsiLeft xs i <|> xsiRight xs i

max_stars_repo_name: eligottlieb/Shannon

max_issues_repo_name: eligottlieb/Shannon

max_forks_repo_name: eligottlieb/Shannon


text:
{- |
This module provides normalized versions of the transforms in @fftw@.

All of the transforms are normalized so that

 - Each transform is unitary, i.e., preserves the inner product and the sum-of-squares norm of its input.

 - Each backwards transform is the inverse of the corresponding forwards transform.

(Both conditions only hold approximately, due to floating point precision.)

For more information on the underlying transforms, see
<http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html>.
-}
module Numeric.FFT.Vector.Unitary(
                -- * Creating and executing 'Plan's
                run,
                plan,
                execute,
                -- * Complex-to-complex transforms
                dft,
                idft,
                -- * Real-to-complex transforms
                dftR2C,
                dftC2R,
                -- * Discrete cosine transforms
                -- $dct_size
                dct2,
                idct2,
                dct4,
                ) where

import Numeric.FFT.Vector.Base
import qualified Numeric.FFT.Vector.Unnormalized as U
import Data.Complex
import qualified Data.Vector.Storable.Mutable as MS
import Control.Monad.Primitive(RealWorld)

-- | A discrete Fourier transform. The output and input sizes are the same (@n@).
--
-- @y_k = (1\/sqrt n) sum_(j=0)^(n-1) x_j e^(-2pi i j k\/n)@
dft :: Transform (Complex Double) (Complex Double)
dft = U.dft {normalization = \n -> constMultOutput $ 1 / sqrt (toEnum n)}

-- | An inverse discrete Fourier transform.  The output and input sizes are the same (@n@).
--
-- @y_k = (1\/sqrt n) sum_(j=0)^(n-1) x_j e^(2pi i j k\/n)@
idft :: Transform (Complex Double) (Complex Double)
idft = U.idft {normalization = \n -> constMultOutput $ 1 / sqrt (toEnum n)}

-- | A forward discrete Fourier transform with real data.  If the input size is @n@,
-- the output size will be @n \`div\` 2 + 1@.
dftR2C :: Transform Double (Complex Double)
dftR2C = U.dftR2C {normalization = \n -> modifyOutput $
                    complexR2CScaling (sqrt 2) n
        }

-- | A normalized backward discrete Fourier transform which is the left inverse of
-- 'U.dftR2C'.  (Specifically, @run dftC2R . run dftR2C == id@.)
--
-- This 'Transform' behaves differently than the others:
--
--  - Calling @plan dftC2R n@ creates a 'Plan' whose /output/ size is @n@, and whose
--    /input/ size is @n \`div\` 2 + 1@.
--
--  - If @length v == n@, then @length (run dftC2R v) == 2*(n-1)@.
--
dftC2R :: Transform (Complex Double) Double
dftC2R = U.dftC2R {normalization = \n -> modifyInput $
                    complexR2CScaling (sqrt 0.5) n
        }

complexR2CScaling :: Double -> Int -> MS.MVector RealWorld (Complex Double) -> IO ()
complexR2CScaling !t !n !a = do
    let !s1 = sqrt (1/toEnum n)
    let !s2 = t * s1
    let len = MS.length a
    -- Justification for the use of unsafeModify:
    -- The output size is 2n+1; so if n>0 then the output size is >=1;
    -- and if n even then the output size is >=3.
    unsafeModify a 0 $ scaleByD s1
    if odd n
        then multC s2 (MS.unsafeSlice 1 (len-1) a)
        else do
            unsafeModify a (len-1) $ scaleByD s1
            multC s2 (MS.unsafeSlice 1 (len-2) a)

-- $dct_size
-- Some normalized real-even (DCT).  The input and output sizes
-- are the same (@n@).

-- | A type-4 discrete cosine transform.  It is its own inverse.
--
-- @y_k = (1\/sqrt n) sum_(j=0)^(n-1) x_j cos(pi(j+1\/2)(k+1\/2)\/n)@
dct4 :: Transform Double Double
dct4 = U.dct4 {normalization = \n -> constMultOutput $ 1 / sqrt (2 * toEnum n)}

-- | A type-2 discrete cosine transform.  Its inverse is 'dct3'.
--
-- @y_k = w(k) sum_(j=0)^(n-1) x_j cos(pi(j+1\/2)k\/n);@
-- where
-- @w(0)=1\/sqrt n@, and @w(k)=sqrt(2\/n)@ for @k>0@.
dct2 :: Transform Double Double
dct2 = U.dct2 {normalization = \n -> modifyOutput $ \a -> do
    let n' = toEnum n
    let !s1 = sqrt $ 1 / (4*n')
    let !s2 = sqrt $ 1 / (2*n')
    unsafeModify a 0 (*s1)
    multC s2 (MS.unsafeSlice 1 (MS.length a-1) a)
    }

-- | A type-3 discrete cosine transform which is the inverse of 'dct2'.
--
-- @y_k = (-1)^k w(n-1) x_(n-1) + 2 sum_(j=0)^(n-2) w(j) x_j sin(pi(j+1)(k+1\/2)/n);@
-- where
-- @w(0)=1\/sqrt(n)@, and @w(k)=1/sqrt(2n)@ for @k>0@.
idct2 :: Transform Double Double
idct2 = U.dct3 {normalization = \n -> modifyInput $ \a -> do
    let n' = toEnum n
    let !s1 = sqrt $ 1 / n'
    let !s2 = sqrt $ 1 / (2*n')
    unsafeModify a 0 (*s1)
    multC s2 (MS.unsafeSlice 1 (MS.length a-1) a)
    }

max_stars_repo_name: TravisWhitaker/vector-fftw

max_issues_repo_name: TravisWhitaker/vector-fftw

max_forks_repo_name: TravisWhitaker/vector-fftw


text:
module GradientDescent where
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Util

--  Performs gradient descent to learn theta
--  theta = GRADIENTDESCENT(x, y, theta, alpha, num_iters) updates theta by
--  taking num_iters gradient steps with learning rate alpha
gradientDescent :: (Matrix Double -> Matrix Double) -> Matrix Double -> Double -> Int -> Matrix Double
gradientDescent gradFunc initialTheta alpha num_iters = iterate initialTheta num_iters
            where
                iterate theta 0 = theta
                iterate theta iters = iterate (updateTheta theta) (iters-1)
                updateTheta t = t - (scalar alpha) * (gradFunc t)

max_stars_repo_name: thade/haskellml

max_issues_repo_name: thade/haskellml

max_forks_repo_name: thade/haskellml


text:
module Scene.Object.SimpleSphere where

import Scene.Base
import Math.Color
import Math.Ray
import Numeric.LinearAlgebra

-- Sphere
sphere :: Double -> Vec3 -> Material -> Scene
sphere r p mat = [sphere' r p mat]

sphere' :: Double -> Vec3 -> Material -> Object
sphere' r p = Object (sphere_intercept r p)

sphere_intercept :: Double -> Vec3 -> Ray -> Maybe (Double, Ray)
sphere_intercept r p ray@(Ray o v) =
    let
        o' = o-p
        a = v <.> v
        b = v <.> o'
        c = o' <.> o' - r**2
        discriminant = b**2 - a*c
    in if discriminant > 0
        then
            let
                x1 = (-b + sqrt(discriminant))/a
                x2 = (-b - sqrt(discriminant))/a
                x = case (x1 > 0, x2 > 0) of
                    (True,True) -> Just (min x1 x2)
                    (False,True) -> Just x2
                    (True, False) -> Just x1
                    (False, False) -> Nothing

            in (\lambda -> (lambda, Ray (parametric ray lambda) (parametric ray lambda - p))) <$> x
        else Nothing

max_stars_repo_name: danrocag/raytracer

max_issues_repo_name: danrocag/raytracer

max_forks_repo_name: danrocag/raytracer


text:
import Mnist
import Runner
import Backprop

import Numeric.LinearAlgebra
import Data.List
import System.Random

learningRate = 0.007

smallRandoms :: Int -> [Double]
smallRandoms seed = map (/100) (randoms (mkStdGen seed))

randomWeightMatrix :: Int -> Int -> Int -> Matrix Double
randomWeightMatrix numInputs numOutputs seed = (numOutputs><numInputs) weights
    where weights = take (numOutputs*numInputs) (smallRandoms seed)

zeroWeightMatrix :: Int -> Int -> Matrix Double
zeroWeightMatrix numInputs numOutputs = (numOutputs><numInputs) weights
    where weights = repeat 0

main :: IO ()
main = do
  let w1 = randomWeightMatrix (28*28 + 1) 20 7
  let w2 = randomWeightMatrix 20 10 42
  let initialNet = buildBackpropNet learningRate [w1, w2] tanhAS
  trainingData2 <- readTrainingData
  let trainingData = take 20000 trainingData2
  putStrLn $ "Training with " ++ show (length trainingData) ++ " images"
  let finalNet = trainWithAllPatterns initialNet trainingData
  testData2 <- readTestData
  let testData = take 1000 testData2
  putStrLn $ "Testing with " ++ show (length testData) ++ " images"
  let results = evalAllPatterns finalNet testData
  let score = fromIntegral (sum results)
  let count = fromIntegral (length testData)
  let percentage = 100.0 * score / count
  putStrLn $ "I got " ++ show percentage ++ "% correct"

max_stars_repo_name: mhwombat/backprop-example

max_issues_repo_name: mhwombat/backprop-example

max_forks_repo_name: mhwombat/backprop-example


text:
module Lambda.Types where

-- Global Imports
import Data.Ratio
import Data.Complex
import Control.Monad.Except
import Text.ParserCombinators.Parsec.Error
import Data.IORef
import System.IO
import qualified Data.Map.Strict as Map

data LispVal
  = Atom String
  | List [LispVal]
  | DottedList [LispVal] LispVal
  | Vector [LispVal]
  | Number Integer
  | Float Double
  | Ratio Rational
  | Complex (Complex Double)
  | String String
  | Bool Bool
  | Character Char
  | IOFunc ([LispVal] -> IOThrowsError LispVal)
  | Port Handle
  | PrimitiveFunc ([LispVal] -> ThrowsError LispVal)
  | Func { params  :: [String]
         , vararg  :: (Maybe String)
         , body    :: [LispVal]
         , closure :: Env
         }

instance Show LispVal where show = showVal
-- This feels like a hack...not sure why it is working or if it truly is
instance Eq LispVal where (==) = (==)

showVal :: LispVal -> String
showVal (String contents)      = "\"" ++ contents ++ "\""
showVal (Atom name)            = name
showVal (Number contents)      = show contents
showVal (Bool True)            = "#t"
showVal (Bool False)           = "#f"
showVal (List contents)        = "(" ++ unwordsList contents ++ ")"
showVal (DottedList head tail) = "(" ++ unwordsList head ++ " . " ++ showVal tail ++ ")"
showVal (Vector contents)      = "(" ++ unwordsList contents ++ ")"
showVal (IOFunc _)             = "<IO primitive>"
showVal (Port _)               = "<IO port>"
showVal (PrimitiveFunc _)      = "<primitive>"
showVal (Func {params = args ,
vararg = varargs             ,
body = body                  ,
               closure = env}) = "(lambda (" ++ unwords (map show args) ++
                 (case varargs of
                    Nothing -> ""
                    Just arg -> " . " ++ arg) ++ ") ...)"

unwordsList :: [LispVal] -> String
unwordsList = unwords . map showVal

-- LispError Data Type and associated functions
data LispError = NumArgs Integer [LispVal]
               | TypeMismatch String LispVal
               | Parser ParseError
               | BadSpecialForm String LispVal
               | NotFunction String String
               | UnboundVar String String
               | Default String

showError :: LispError -> String
showError (UnboundVar msg var)          = msg ++ ": " ++ var
showError (BadSpecialForm msg form)     = msg ++ ": " ++ show form
showError (NotFunction msg func)        = msg ++ ": " ++ show func
showError (NumArgs expected found)      = "Expected " ++ show expected ++ " args; found values " ++ unwordsList found
showError (TypeMismatch expected found) = "Invalid type.  Expected " ++ expected ++ ", found " ++ show found
showError (Parser parseErr)             = "Parse error at " ++ show parseErr
showError (Default err)                 = "Unknown error: " ++ show err

instance Show LispError where show = showError

type ThrowsError = Either LispError

-- Helper Functions for Parsing
trapError action = catchError action (return . show)

extractValue :: ThrowsError a -> a
extractValue (Right val) = val

-- Environment Errors
type IOThrowsError = ExceptT LispError IO

liftThrows :: ThrowsError a -> IOThrowsError a
-- Note that this type signature is not general...
-- try (MonadError m a) => Either e a -> m a
liftThrows (Left err)  = throwError err
liftThrows (Right val) = return val

runIOThrows :: IOThrowsError String -> IO String
runIOThrows action = runExceptT (trapError action) >>= return . extractValue

-- Environment Definition
type Env = IORef (Map.Map String (IORef LispVal))

nullEnv :: IO Env
nullEnv = newIORef Map.empty

isBound :: Env -> String -> IO Bool
isBound e s = readIORef e >>= return . Map.member s

getVar :: Env -> String -> IOThrowsError LispVal
getVar e s = do
  env <- liftIO $ readIORef e
  case Map.lookup s env of
    Just v  -> liftIO $ readIORef v
    Nothing -> throwError $ UnboundVar "Trying to get an unbound variable..." s

setVar :: Env -> String -> LispVal -> IOThrowsError LispVal
setVar e s l = do
  env <- liftIO $ readIORef e
  case Map.lookup s env of
    Just v  -> liftIO $ writeIORef v l
    Nothing -> throwError $ UnboundVar "Trying to update an unbound variable.." s
  return l

defineVar :: Env -> String -> LispVal -> IOThrowsError LispVal
defineVar e s l = do
  defined <- liftIO $ isBound e s
  if defined then setVar e s l
  else liftIO $ do
    newVar <- newIORef l
    modifyIORef' e $ Map.insert s newVar
    return l

bindVars :: Env -> [(String, LispVal)] -> IO Env
bindVars e ls = do
  updates <- mapM mkRef ls >>= return . Map.fromList
  modifyIORef' e (Map.union updates)
  return e
    where
      mkRef :: (String, LispVal) -> IO (String, IORef LispVal)
      mkRef (s, l) = newIORef l >>= return . (,) s

showEnv :: Env -> IO String
showEnv e = do
  env <- readIORef e
  vars <- return $ Map.toList env
  (sequence $ map printOne vars) >>= return . unlines
    where
      printOne :: (String, IORef LispVal) -> IO String
      printOne (s, l) = do
        lispval <- readIORef l
        return $ s ++ " : " ++ show lispval

max_stars_repo_name: tdfirth/lambda

max_issues_repo_name: tdfirth/lambda

max_forks_repo_name: tdfirth/lambda


text:
{-# LANGUAGE PatternGuards #-}
-- |
-- Module    : Statistics.Matrix
-- Copyright : 2011 Aleksey Khudyakov, 2014 Bryan O'Sullivan
-- License   : BSD3
--
-- Basic matrix operations.
--
-- There isn't a widely used matrix package for Haskell yet, so
-- we implement the necessary minimum here.

module Statistics.Matrix
    ( -- * Data types
      Matrix(..)
    , Vector
      -- * Conversion from/to lists/vectors
    , fromVector
    , fromList
    , fromRowLists
    , fromRows
    , fromColumns
    , toVector
    , toList
    , toRows
    , toColumns
    , toRowLists
      -- * Other
    , generate
    , generateSym
    , ident
    , diag
    , dimension
    , center
    , multiply
    , multiplyV
    , transpose
    , power
    , norm
    , column
    , row
    , map
    , for
    , unsafeIndex
    , hasNaN
    , bounds
    , unsafeBounds
    ) where

import Prelude hiding (exponent, map)
import Control.Applicative ((<$>))
import Control.Monad.ST
import qualified Data.Vector.Unboxed as U
import           Data.Vector.Unboxed   ((!))
import qualified Data.Vector.Unboxed.Mutable as UM

import Numeric.Sum         (sumVector,kbn)
import Statistics.Matrix.Function
import Statistics.Matrix.Types
import Statistics.Matrix.Mutable  (unsafeNew,unsafeWrite,unsafeFreeze)

----------------------------------------------------------------
-- Conversion to/from vectors/lists
----------------------------------------------------------------

-- | Convert from a row-major list.
fromList :: Int                 -- ^ Number of rows.
         -> Int                 -- ^ Number of columns.
         -> [Double]            -- ^ Flat list of values, in row-major order.
         -> Matrix
fromList r c = fromVector r c . U.fromList

-- | create a matrix from a list of lists, as rows
fromRowLists :: [[Double]] -> Matrix
fromRowLists = fromRows . fmap U.fromList

-- | Convert from a row-major vector.
fromVector :: Int               -- ^ Number of rows.
           -> Int               -- ^ Number of columns.
           -> U.Vector Double   -- ^ Flat list of values, in row-major order.
           -> Matrix
fromVector r c v
  | r*c /= len = error "input size mismatch"
  | otherwise  = Matrix r c v
  where len    = U.length v

-- | create a matrix from a list of vectors, as rows
fromRows :: [Vector] -> Matrix
fromRows xs
  | [] <- xs        = error "Statistics.Matrix.fromRows: empty list of rows!"
  | any (/=nCol) ns = error "Statistics.Matrix.fromRows: row sizes do not match"
  | nCol == 0       = error "Statistics.Matrix.fromRows: zero columns in matrix"
  | otherwise       = fromVector nRow nCol (U.concat xs)
  where
    nCol:ns = U.length <$> xs
    nRow    = length xs

-- | create a matrix from a list of vectors, as columns
fromColumns :: [Vector] -> Matrix
fromColumns = transpose . fromRows

-- | Convert to a row-major flat vector.
toVector :: Matrix -> U.Vector Double
toVector (Matrix _ _ v) = v

-- | Convert to a row-major flat list.
toList :: Matrix -> [Double]
toList = U.toList . toVector

-- | Convert to a list of lists, as rows
toRowLists :: Matrix -> [[Double]]
toRowLists (Matrix _ nCol v)
  = chunks $ U.toList v
  where
    chunks [] = []
    chunks xs = case splitAt nCol xs of
      (rowE,rest) -> rowE : chunks rest

-- | Convert to a list of vectors, as rows
toRows :: Matrix -> [Vector]
toRows (Matrix _ nCol v) = chunks v
  where
    chunks xs
      | U.null xs = []
      | otherwise = case U.splitAt nCol xs of
          (rowE,rest) -> rowE : chunks rest

-- | Convert to a list of vectors, as columns
toColumns :: Matrix -> [Vector]
toColumns = toRows . transpose

----------------------------------------------------------------
-- Other
----------------------------------------------------------------

-- | Generate matrix using function
generate :: Int                 -- ^ Number of rows
         -> Int                 -- ^ Number of columns
         -> (Int -> Int -> Double)
            -- ^ Function which takes /row/ and /column/ as argument.
         -> Matrix
generate nRow nCol f
  = Matrix nRow nCol $ U.generate (nRow*nCol) $ \i ->
      let (r,c) = i `quotRem` nCol in f r c

-- | Generate symmetric square matrix using function
generateSym
  :: Int                 -- ^ Number of rows and columns
  -> (Int -> Int -> Double)
     -- ^ Function which takes /row/ and /column/ as argument. It must
     --   be symmetric in arguments: @f i j == f j i@
  -> Matrix
generateSym n f = runST $ do
  m <- unsafeNew n n
  for 0 n $ \r -> do
    unsafeWrite m r r (f r r)
    for (r+1) n $ \c -> do
      let x = f r c
      unsafeWrite m r c x
      unsafeWrite m c r x
  unsafeFreeze m

-- | Create the square identity matrix with given dimensions.
ident :: Int -> Matrix
ident n = diag $ U.replicate n 1.0

-- | Create a square matrix with given diagonal, other entries default to 0
diag :: Vector -> Matrix
diag v
  = Matrix n n $ U.create $ do
      arr <- UM.replicate (n*n) 0
      for 0 n $ \i ->
        UM.unsafeWrite arr (i*n + i) (v ! i)
      return arr
  where
    n = U.length v

-- | Return the dimensions of this matrix, as a (row,column) pair.
dimension :: Matrix -> (Int, Int)
dimension (Matrix r c _) = (r, c)

-- | Matrix-matrix multiplication. Matrices must be of compatible
-- sizes (/note: not checked/).
multiply :: Matrix -> Matrix -> Matrix
multiply m1@(Matrix r1 _ _) m2@(Matrix _ c2 _) =
  Matrix r1 c2 $ U.generate (r1*c2) go
  where
    go t = sumVector kbn $ U.zipWith (*) (row m1 i) (column m2 j)
      where (i,j) = t `quotRem` c2

-- | Matrix-vector multiplication.
multiplyV :: Matrix -> Vector -> Vector
multiplyV m v
  | cols m == c = U.generate (rows m) (sumVector kbn . U.zipWith (*) v . row m)
  | otherwise   = error $ "matrix/vector unconformable " ++ show (cols m,c)
  where c = U.length v

-- | Raise matrix to /n/th power. Power must be positive
-- (/note: not checked).
power :: Matrix -> Int -> Matrix
power mat 1 = mat
power mat n = res
  where
    mat2 = power mat (n `quot` 2)
    pow  = multiply mat2 mat2
    res | odd n     = multiply pow mat
        | otherwise = pow

-- | Element in the center of matrix (not corrected for exponent).
center :: Matrix -> Double
center mat@(Matrix r c _) =
    unsafeBounds U.unsafeIndex mat (r `quot` 2) (c `quot` 2)

-- | Calculate the Euclidean norm of a vector.
norm :: Vector -> Double
norm = sqrt . sumVector kbn . U.map square

-- | Return the given column.
column :: Matrix -> Int -> Vector
column (Matrix r c v) j= U.generate r (\i -> v `U.unsafeIndex` (j + i * c))
{-# INLINE column #-}

-- | Return the given row.
row :: Matrix -> Int -> Vector
row (Matrix _ c v) i = U.slice (c*i) c v

unsafeIndex :: Matrix
            -> Int              -- ^ Row.
            -> Int              -- ^ Column.
            -> Double
unsafeIndex = unsafeBounds U.unsafeIndex

-- | Apply function to every element of matrix
map :: (Double -> Double) -> Matrix -> Matrix
map f (Matrix r c v) = Matrix r c (U.map f v)

-- | Indicate whether any element of the matrix is @NaN@.
hasNaN :: Matrix -> Bool
hasNaN = U.any isNaN . toVector

-- | Given row and column numbers, calculate the offset into the flat
-- row-major vector.
bounds :: (Vector -> Int -> r) -> Matrix -> Int -> Int -> r
bounds k (Matrix rs cs v) r c
  | r < 0 || r >= rs = error "row out of bounds"
  | c < 0 || c >= cs = error "column out of bounds"
  | otherwise        = k v $! r * cs + c
{-# INLINE bounds #-}

-- | Given row and column numbers, calculate the offset into the flat
-- row-major vector, without checking.
unsafeBounds :: (Vector -> Int -> r) -> Matrix -> Int -> Int -> r
unsafeBounds k (Matrix _ cs v) r c = k v $! r * cs + c
{-# INLINE unsafeBounds #-}

transpose :: Matrix -> Matrix
transpose m@(Matrix r0 c0 _) = Matrix c0 r0 . U.generate (r0*c0) $ \i ->
  let (r,c) = i `quotRem` r0
  in unsafeIndex m c r

max_stars_repo_name: arvindd/dh-core

max_issues_repo_name: arvindd/dh-core

max_forks_repo_name: arvindd/dh-core


text:
module Approx ( bezier2biarc
              ) where

import qualified CubicBezier as B
import qualified BiArc as BA
import qualified Line as L

import Data.Bool (bool)
import Linear
import Data.Complex

import Types

-- Approximate a bezier curve with biarcs (Left) and line segments (Right)
bezier2biarc :: B.CubicBezier
             -> Double
             -> [Either BA.BiArc (V2 Double)]
bezier2biarc mbezier resolution
    -- Edge case: all points on the same line -> it is a line
    | (L.isOnLine (L.fromPoints (B._p2 mbezier) (B._p1 mbezier)) (B._c1 mbezier)) &&
      (L.isOnLine (L.fromPoints (B._p2 mbezier) (B._p1 mbezier)) (B._c2 mbezier))
        = [Right (B._p2 mbezier)]
    -- Edge case: p1 == c1, don't split
    | (B._p1 mbezier) == (B._c1 mbezier)
        = approxOne mbezier
    -- Edge case: p2 == c2, don't split
    | (B._p2 mbezier) == (B._c2 mbezier)
        = approxOne mbezier
    -- Split by the inflexion points (if any)
    | otherwise
        = byInflection (B.realInflectionPoint i1) (B.realInflectionPoint i2)
    where
        (i1, i2) = B.inflectionPoints mbezier

        order a b | b < a = (b, a)
                  | otherwise = (a, b)

        byInflection True False = approxOne b1 ++ approxOne b2
            where
                (b1, b2) = B.bezierSplitAt mbezier (realPart i1)

        byInflection False True = approxOne b1 ++ approxOne b2
            where
                (b1, b2) = B.bezierSplitAt mbezier (realPart i2)

        byInflection True True = approxOne b1 ++ approxOne b2 ++ approxOne b3
            where
                (it1, it2') = order (realPart i1) (realPart i2)

                -- Make the first split and save the first new curve. The second one has to be splitted again
                -- at the recalculated t2 (it is on a new curve)
                it2 = (1 - it1) * it2'

                (b1, toSplit) = B.bezierSplitAt mbezier it1
                (b2, b3) = B.bezierSplitAt toSplit it2

        byInflection False False = approxOne mbezier

        -- TODO: make it tail recursive
        approxOne :: B.CubicBezier -> [Either BA.BiArc (V2 Double)]
        approxOne bezier
            -- Approximate bezier length. if smaller than resolution, do not approximate
            | (distance (B._p1 bezier) (B._c1 bezier)) +
              (distance (B._c1 bezier) (B._c2 bezier)) +
              (distance (B._c2 bezier) (B._p2 bezier)) < resolution
                = [Right (B._p2 bezier)]
            -- Edge case: start- and endpoints are the same
            | (B._p1 bezier) == (B._p2 bezier)
                = splitAndRecur 0.5
            -- Edge case: control lines are parallel
            | (L._m t1) == (L._m t2) || (isNaN (L._m t1) && isNaN (L._m t2))
                = splitAndRecur 0.5
            -- Approximation is not close enough yet, refine
            | BA.isStable biarc && maxDistance > resolution
                = splitAndRecur maxDistanceAt
            -- Desired case: approximation is stable and close enough
            | BA.isStable biarc
                = [Left biarc]
            -- Unstable approximation: split the bezier into half, basically switching to
            -- linear approximation mode
            | otherwise
                = splitAndRecur 0.5

            where
                -- Edge case: P1==C1 or P2==C2
                -- there is no derivative at P1 or P2, use the other control point
                c1 = bool (B._c1 bezier) (B._c2 bezier) ((B._p1 bezier) == (B._c1 bezier))
                c2 = bool (B._c2 bezier) (B._c1 bezier) ((B._p2 bezier) == (B._c2 bezier))

                -- V: Intersection point of tangent lines
                t1 = L.fromPoints (B._p1 bezier) c1
                t2 = L.fromPoints (B._p2 bezier) c2
                v = L.intersection t1 t2

                -- G: incenter point of the triangle (P1, V, P2)
                dP2V = distance (B._p2 bezier) v
                dP1V = distance (B._p1 bezier) v
                dP1P2 = distance (B._p1 bezier) (B._p2 bezier)
                g = (dP2V *^ B._p1 bezier + dP1V *^ B._p2 bezier + dP1P2 *^ v) ^/ (dP2V + dP1V + dP1P2)

                -- Calculate the BiArc
                biarc = BA.create (B._p1 bezier) (B._p1 bezier - c1) (B._p2 bezier) (B._p2 bezier - c2) g

                -- Calculate the error
                -- TODO: we only calculate the distance at 8 points (first and last skipped as
                --       they should be precise), seems a resonable approximation as for now
                parameterStep = 1 / 10

                (maxDistance, maxDistanceAt) = maxDistance' 0 0 parameterStep

                maxDistance' m mt t
                    | t < 1
                        = if' (d > m) (maxDistance' d t nt) (maxDistance' m mt nt)
                    | otherwise
                        = (m, mt)
                    where
                        d = distance (BA.pointAt biarc t) (B.pointAt bezier t)
                        nt = t + parameterStep

                splitAndRecur t = let (b1, b2) = B.bezierSplitAt bezier t
                                   in approxOne b1 ++ approxOne b2


text:
module FastFourier (fastFourier, inverseFastFourier, getNewIdxOrder) where

import Data.Complex (Complex, cis, magnitude)
import Data.List (head, tail, foldl, zipWith, length, reverse)
import Data.List.Split (chunksOf)
import Data.Char (intToDigit, digitToInt)
import Numeric (showIntAtBase)

import ComplexUtils (toComplexM)

splitByOddAndEvenIdx :: [(Complex Double)] -> ([(Complex Double)], [(Complex Double)])
splitByOddAndEvenIdx items =
    let
        chunks = chunksOf 2 items
        even = map head chunks
        odd = map last chunks
    in
        (even, odd)

splitInHalf :: [(Complex Double)] -> ([(Complex Double)], [(Complex Double)])
splitInHalf items =
    let l = (length items) `div` 2
    in (take l items, drop l items)

getAngle :: Int -> Int -> Int -> (Complex Double)
getAngle op m k =
    if k == 0
        then cis 0
        else cis (2.0 * pi * (fromIntegral k) * (fromIntegral op) / (fromIntegral m))

getAngles :: Int -> Int -> [Int] -> [(Complex Double)]
getAngles op m ks = map (getAngle op m) ks

fastFourier' :: Int -> [(Complex Double)] -> [(Complex Double)]
fastFourier' op fs =
    case fs of
        [f] -> [f]
        fs ->
            let
                m = length fs
                k_max = (m `div` 2) - 1
                ks = [0..k_max]

--                 (even, odd) = splitByOddAndEvenIdx(fs)
                (even, odd) = splitInHalf(fs)

                a_s = fastFourier' op even
                b_s = fastFourier' op odd

                angles = getAngles op m ks

                w_b_s = zipWith (*) b_s angles
                c_head = zipWith (+) a_s w_b_s
                c_tail = zipWith (-) a_s w_b_s
            in
                c_head ++ c_tail

addZeros :: Num a => Int -> [a] -> [a]
addZeros maxSize arr =
        let diff = maxSize - (length arr)
        in arr ++ replicate diff 0

padWithZeros :: [Double] -> [Double]
padWithZeros fs =
    let
        len = length fs
        pow = ceiling $ logBase 2 $ fromIntegral len
        a = 2^pow
    in
        addZeros a fs

toInt :: String -> Int
toInt s = foldl (\acc x -> acc * 2 + digitToInt x) 0 s

toBinaryStr :: Int -> String
toBinaryStr i = showIntAtBase 2 intToDigit i ""

addStrZeros :: Int -> String -> String
addStrZeros maxSize arr =
        let diff = maxSize - (length arr)
        in arr ++ replicate diff '0'

binaryInverse :: Int -> Int -> Int
binaryInverse maxDigitCount i =
    toInt $ addStrZeros maxDigitCount $ reverse $ toBinaryStr i

getNewIdxOrder :: Int -> [Int]
getNewIdxOrder len =
    let
        n = len -1
        maxDigitCnt = length $ toBinaryStr n
    in
        map (binaryInverse maxDigitCnt) [0..n]

reorder :: [a] -> [a]
reorder fs =
    map (\x -> fs!!x) $ getNewIdxOrder $ length fs

fastFourier :: [Double] -> [(Complex Double)]
fastFourier fs = fastFourier' (-1) $ toComplexM $ reorder $ padWithZeros fs

prepareForInverse :: [(Complex Double)] -> [(Complex Double)]
prepareForInverse cs =
    let c0 = head cs
    in [c0] ++ (reverse $ tail cs)

inverseFastFourier :: [(Complex Double)] -> [Double]
inverseFastFourier cs =
    let
        n = length cs
        nd = fromIntegral n
    in
        map (/nd) $ map magnitude $ fastFourier' (1) $ reorder cs

max_stars_repo_name: vined/fourier

max_issues_repo_name: vined/fourier

max_forks_repo_name: vined/fourier