text:
{-# LANGUAGE CPP                #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE DeriveFoldable     #-}
{-# LANGUAGE DeriveFunctor      #-}
{-# LANGUAGE DeriveGeneric      #-}
{-# LANGUAGE FlexibleContexts   #-}
{-# LANGUAGE ViewPatterns       #-}
-- |
-- Module    : Statistics.Quantile
-- Copyright : (c) 2009 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Functions for approximating quantiles, i.e. points taken at regular
-- intervals from the cumulative distribution function of a random
-- variable.
--
-- The number of quantiles is described below by the variable /q/, so
-- with /q/=4, a 4-quantile (also known as a /quartile/) has 4
-- intervals, and contains 5 points.  The parameter /k/ describes the
-- desired point, where 0 ≤ /k/ ≤ /q/.

module Statistics.Quantile
    (
    -- * Quantile estimation functions
    -- $cont_quantiles
      ContParam(..)
    , Default(..)
    , quantile
    , quantiles
    , quantilesVec
    -- ** Parameters for the continuous sample method
    , cadpw
    , hazen
    , spss
    , s
    , medianUnbiased
    , normalUnbiased
    -- * Other algorithms
    , weightedAvg
    -- * Median & other specializations
    , median
    , mad
    , midspread
    -- * Deprecated
    , continuousBy
    -- * References
    -- $references
    ) where

import           Data.Binary            (Binary)
import           Data.Aeson             (ToJSON,FromJSON)
import           Data.Data              (Data,Typeable)
import           Data.Default.Class
import           Data.Functor
import qualified Data.Foldable        as F
import           Data.Vector.Generic ((!))
import qualified Data.Vector          as V
import qualified Data.Vector.Generic  as G
import qualified Data.Vector.Unboxed  as U
import qualified Data.Vector.Storable as S
import GHC.Generics (Generic)

import Statistics.Function (partialSort)

----------------------------------------------------------------
-- Quantile estimation
----------------------------------------------------------------

-- | O(/n/·log /n/). Estimate the /k/th /q/-quantile of a sample,
-- using the weighted average method. Up to rounding errors it's same
-- as @quantile s@.
--
-- The following properties should hold otherwise an error will be thrown.
--
--   * the length of the input is greater than @0@
--
--   * the input does not contain @NaN@
--
--   * k ≥ 0 and k ≤ q
weightedAvg :: G.Vector v Double =>
               Int        -- ^ /k/, the desired quantile.
            -> Int        -- ^ /q/, the number of quantiles.
            -> v Double   -- ^ /x/, the sample data.
            -> Double
weightedAvg k q x
  | G.any isNaN x   = modErr "weightedAvg" "Sample contains NaNs"
  | n == 0          = modErr "weightedAvg" "Sample is empty"
  | n == 1          = G.head x
  | q < 2           = modErr "weightedAvg" "At least 2 quantiles is needed"
  | k == q          = G.maximum x
  | k >= 0 || k < q = xj + g * (xj1 - xj)
  | otherwise       = modErr "weightedAvg" "Wrong quantile number"
  where
    j   = floor idx
    idx = fromIntegral (n - 1) * fromIntegral k / fromIntegral q
    g   = idx - fromIntegral j
    xj  = sx ! j
    xj1 = sx ! (j+1)
    sx  = partialSort (j+2) x
    n   = G.length x
{-# SPECIALIZE weightedAvg :: Int -> Int -> U.Vector Double -> Double #-}
{-# SPECIALIZE weightedAvg :: Int -> Int -> V.Vector Double -> Double #-}
{-# SPECIALIZE weightedAvg :: Int -> Int -> S.Vector Double -> Double #-}

----------------------------------------------------------------
-- Quantiles continuous algorithm
----------------------------------------------------------------

-- $cont_quantiles
--
-- Below is family of functions which use same algorithm for estimation
-- of sample quantiles. It approximates empirical CDF as continuous
-- piecewise function which interpolates linearly between points
-- \((X_k,p_k)\) where \(X_k\) is k-th order statistics (k-th smallest
-- element) and \(p_k\) is probability corresponding to
-- it. 'ContParam' determines how \(p_k\) is chosen. For more detailed
-- explanation see [Hyndman1996].
--
-- This is the method used by most statistical software, such as R,
-- Mathematica, SPSS, and S.

-- | Parameters /α/ and /β/ to the 'continuousBy' function. Exact
--   meaning of parameters is described in [Hyndman1996] in section
--   \"Piecewise linear functions\"
data ContParam = ContParam {-# UNPACK #-} !Double {-# UNPACK #-} !Double
  deriving (Show,Eq,Ord,Data,Typeable,Generic)

-- | We use 's' as default value which is same as R's default.
instance Default ContParam where
  def = s

instance Binary   ContParam
instance ToJSON   ContParam
instance FromJSON ContParam

-- | O(/n/·log /n/). Estimate the /k/th /q/-quantile of a sample /x/,
--   using the continuous sample method with the given parameters.
--
--   The following properties should hold, otherwise an error will be thrown.
--
--     * input sample must be nonempty
--
--     * the input does not contain @NaN@
--
--     * 0 ≤ k ≤ q
quantile :: G.Vector v Double
         => ContParam  -- ^ Parameters /α/ and /β/.
         -> Int        -- ^ /k/, the desired quantile.
         -> Int        -- ^ /q/, the number of quantiles.
         -> v Double   -- ^ /x/, the sample data.
         -> Double
quantile param q nQ xs
  | nQ < 2         = modErr "continuousBy" "At least 2 quantiles is needed"
  | badQ nQ q      = modErr "continuousBy" "Wrong quantile number"
  | G.any isNaN xs = modErr "continuousBy" "Sample contains NaNs"
  | otherwise      = estimateQuantile sortedXs pk
  where
    pk       = toPk param n q nQ
    sortedXs = psort xs $ floor pk + 1
    n        = G.length xs
{-# INLINABLE quantile #-}
{-# SPECIALIZE
    quantile :: ContParam -> Int -> Int -> U.Vector Double -> Double #-}
{-# SPECIALIZE
    quantile :: ContParam -> Int -> Int -> V.Vector Double -> Double #-}
{-# SPECIALIZE
    quantile :: ContParam -> Int -> Int -> S.Vector Double -> Double #-}

-- | O(/k·n/·log /n/). Estimate set of the /k/th /q/-quantile of a
--   sample /x/, using the continuous sample method with the given
--   parameters. This is faster than calling quantile repeatedly since
--   sample should be sorted only once
--
--   The following properties should hold, otherwise an error will be thrown.
--
--     * input sample must be nonempty
--
--     * the input does not contain @NaN@
--
--     * for every k in set of quantiles 0 ≤ k ≤ q
quantiles :: (G.Vector v Double, F.Foldable f, Functor f)
  => ContParam
  -> f Int
  -> Int
  -> v Double
  -> f Double
quantiles param qs nQ xs
  | nQ < 2             = modErr "quantiles" "At least 2 quantiles is needed"
  | F.any (badQ nQ) qs = modErr "quantiles" "Wrong quantile number"
  | G.any isNaN xs     = modErr "quantiles" "Sample contains NaNs"
  -- Doesn't matter what we put into empty container
  | fnull qs           = 0 <$ qs
  | otherwise          = fmap (estimateQuantile sortedXs) ks'
  where
    ks'      = fmap (\q -> toPk param n q nQ) qs
    sortedXs = psort xs $ floor (F.maximum ks') + 1
    n        = G.length xs
{-# INLINABLE quantiles #-}
{-# SPECIALIZE quantiles
      :: (Functor f, F.Foldable f) => ContParam -> f Int -> Int -> V.Vector Double -> f Double #-}
{-# SPECIALIZE quantiles
      :: (Functor f, F.Foldable f) => ContParam -> f Int -> Int -> U.Vector Double -> f Double #-}
{-# SPECIALIZE quantiles
      :: (Functor f, F.Foldable f) => ContParam -> f Int -> Int -> S.Vector Double -> f Double #-}

-- COMPAT
fnull :: F.Foldable f => f a -> Bool
#if !MIN_VERSION_base(4,8,0)
fnull = F.foldr (\_ _ -> False) True
#else
fnull = null
#endif

-- | O(/k·n/·log /n/). Same as quantiles but uses 'G.Vector' container
--   instead of 'Foldable' one.
quantilesVec :: (G.Vector v Double, G.Vector v Int)
  => ContParam
  -> v Int
  -> Int
  -> v Double
  -> v Double
quantilesVec param qs nQ xs
  | nQ < 2             = modErr "quantilesVec" "At least 2 quantiles is needed"
  | G.any (badQ nQ) qs = modErr "quantilesVec" "Wrong quantile number"
  | G.any isNaN xs     = modErr "quantilesVec" "Sample contains NaNs"
  | G.null qs          = G.empty
  | otherwise          = G.map (estimateQuantile sortedXs) ks'
  where
    ks'      = G.map (\q -> toPk param n q nQ) qs
    sortedXs = psort xs $ floor (G.maximum ks') + 1
    n        = G.length xs
{-# INLINABLE quantilesVec #-}
{-# SPECIALIZE quantilesVec
      :: ContParam -> V.Vector Int -> Int -> V.Vector Double -> V.Vector Double #-}
{-# SPECIALIZE quantilesVec
      :: ContParam -> U.Vector Int -> Int -> U.Vector Double -> U.Vector Double #-}
{-# SPECIALIZE quantilesVec
      :: ContParam -> S.Vector Int -> Int -> S.Vector Double -> S.Vector Double #-}

-- Returns True if quantile number is out of range
badQ :: Int -> Int -> Bool
badQ nQ q = q < 0 || q > nQ

-- Obtain k from equation for p_k [Hyndman1996] p.363.  Note that
-- equation defines p_k for integer k but we calculate it as real
-- value and will use fractional part for linear interpolation. This
-- is correct since equation is linear.
toPk
  :: ContParam
  -> Int        -- ^ /n/ number of elements
  -> Int        -- ^ /k/, the desired quantile.
  -> Int        -- ^ /q/, the number of quantiles.
  -> Double
toPk (ContParam a b) (fromIntegral -> n) q nQ
  = a + p * (n + 1 - a - b)
  where
    p = fromIntegral q / fromIntegral nQ

-- Estimate quantile for given k (including fractional part)
estimateQuantile :: G.Vector v Double => v Double -> Double -> Double
{-# INLINE estimateQuantile #-}
estimateQuantile sortedXs k'
  = (1-g) * item (k-1) + g * item k
  where
    (k,g) = properFraction k'
    item  = (sortedXs !) . clamp
    --
    clamp = max 0 . min (n - 1)
    n     = G.length sortedXs

psort :: G.Vector v Double => v Double -> Int -> v Double
psort xs k = partialSort (max 0 $ min (G.length xs - 1) k) xs
{-# INLINE psort #-}

-- | California Department of Public Works definition, /α/=0, /β/=1.
-- Gives a linear interpolation of the empirical CDF.  This
-- corresponds to method 4 in R and Mathematica.
cadpw :: ContParam
cadpw = ContParam 0 1

-- | Hazen's definition, /α/=0.5, /β/=0.5.  This is claimed to be
-- popular among hydrologists.  This corresponds to method 5 in R and
-- Mathematica.
hazen :: ContParam
hazen = ContParam 0.5 0.5

-- | Definition used by the SPSS statistics application, with /α/=0,
-- /β/=0 (also known as Weibull's definition).  This corresponds to
-- method 6 in R and Mathematica.
spss :: ContParam
spss = ContParam 0 0

-- | Definition used by the S statistics application, with /α/=1,
-- /β/=1.  The interpolation points divide the sample range into @n-1@
-- intervals.  This corresponds to method 7 in R and Mathematica and
-- is default in R.
s :: ContParam
s = ContParam 1 1

-- | Median unbiased definition, /α/=1\/3, /β/=1\/3. The resulting
-- quantile estimates are approximately median unbiased regardless of
-- the distribution of /x/.  This corresponds to method 8 in R and
-- Mathematica.
medianUnbiased :: ContParam
medianUnbiased = ContParam third third
    where third = 1/3

-- | Normal unbiased definition, /α/=3\/8, /β/=3\/8.  An approximately
-- unbiased estimate if the empirical distribution approximates the
-- normal distribution.  This corresponds to method 9 in R and
-- Mathematica.
normalUnbiased :: ContParam
normalUnbiased = ContParam ta ta
    where ta = 3/8

modErr :: String -> String -> a
modErr f err = error $ "Statistics.Quantile." ++ f ++ ": " ++ err

----------------------------------------------------------------
-- Specializations
----------------------------------------------------------------

-- | O(/n/·log /n/) Estimate median of sample
median :: G.Vector v Double
       => ContParam  -- ^ Parameters /α/ and /β/.
       -> v Double   -- ^ /x/, the sample data.
       -> Double
{-# INLINE median #-}
median p = quantile p 1 2

-- | O(/n/·log /n/). Estimate the range between /q/-quantiles 1 and
-- /q/-1 of a sample /x/, using the continuous sample method with the
-- given parameters.
--
-- For instance, the interquartile range (IQR) can be estimated as
-- follows:
--
-- > midspread medianUnbiased 4 (U.fromList [1,1,2,2,3])
-- > ==> 1.333333
midspread :: G.Vector v Double =>
             ContParam  -- ^ Parameters /α/ and /β/.
          -> Int        -- ^ /q/, the number of quantiles.
          -> v Double   -- ^ /x/, the sample data.
          -> Double
midspread param k x
  | G.any isNaN x = modErr "midspread" "Sample contains NaNs"
  | k <= 0        = modErr "midspread" "Nonpositive number of quantiles"
  | otherwise     = let Pair x1 x2 = quantiles param (Pair 1 (k-1)) k x
                    in  x2 - x1
{-# INLINABLE  midspread #-}
{-# SPECIALIZE midspread :: ContParam -> Int -> U.Vector Double -> Double #-}
{-# SPECIALIZE midspread :: ContParam -> Int -> V.Vector Double -> Double #-}
{-# SPECIALIZE midspread :: ContParam -> Int -> S.Vector Double -> Double #-}

data Pair a = Pair !a !a
  deriving (Functor, F.Foldable)

-- | O(/n/·log /n/). Estimate the median absolute deviation (MAD) of a
--   sample /x/ using 'continuousBy'. It's robust estimate of
--   variability in sample and defined as:
--
--   \[
--   MAD = \operatorname{median}(| X_i - \operatorname{median}(X) |)
--   \]
mad :: G.Vector v Double
    => ContParam  -- ^ Parameters /α/ and /β/.
    -> v Double   -- ^ /x/, the sample data.
    -> Double
mad p xs
  = median p $ G.map (abs . subtract med) xs
  where
    med = median p xs
{-# INLINABLE  mad #-}
{-# SPECIALIZE mad :: ContParam -> U.Vector Double -> Double #-}
{-# SPECIALIZE mad :: ContParam -> V.Vector Double -> Double #-}
{-# SPECIALIZE mad :: ContParam -> S.Vector Double -> Double #-}

----------------------------------------------------------------
-- Deprecated
----------------------------------------------------------------

continuousBy :: G.Vector v Double =>
                ContParam  -- ^ Parameters /α/ and /β/.
             -> Int        -- ^ /k/, the desired quantile.
             -> Int        -- ^ /q/, the number of quantiles.
             -> v Double   -- ^ /x/, the sample data.
             -> Double
continuousBy = quantile
{-# DEPRECATED continuousBy "Use quantile instead" #-}

-- $references
--
-- * Weisstein, E.W. Quantile. /MathWorld/.
--   <http://mathworld.wolfram.com/Quantile.html>
--
-- * [Hyndman1996] Hyndman, R.J.; Fan, Y. (1996) Sample quantiles in statistical
--   packages. /American Statistician/
--   50(4):361–365. <http://www.jstor.org/stable/2684934>


text:
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}

-- | TODO
module SatOpt.Utility where

import           Data.Packed.Matrix            (buildMatrix, liftMatrix2)
import           Data.Packed.Vector            (dim, zipVectorWith)
import           Numeric.LinearAlgebra.Data
import           Numeric.LinearAlgebra.HMatrix

zipMatrixWith ::(Element a, Element b, Element c)
                => (a -> b -> c)
                -> Matrix a
                -> Matrix b
                -> Matrix c
zipMatrixWith f = liftMatrix2 (zipVectorWith f)

sample2D :: (a -> a -> b) -> [a] -> [a] -> [[b]]
sample2D f xs ys = [[f x y | y <- ys] | x <- xs]

sample2DM :: (Element a, Element b, Indexable (Vector a) a)
             => (a -> a -> b)
             -> Vector a
             -> Vector a
             -> Matrix b
sample2DM f vx vy = buildMatrix (dim vx) (dim vy) (\(i, j) -> f (vx ! i) (vy ! j))

square :: Num a => a -> a
square x = x * x

clamp :: Double -> Double
clamp x
  | x > 1     = 1
  | x < 0     = 0
  | otherwise = x

norm :: Floating a => a -> a -> a
norm x y = sqrt $ square x + square y

ext: hs

max_stars_repo_name: taktoa/SatOpt

max_issues_repo_name: taktoa/SatOpt

max_forks_repo_name: taktoa/SatOpt


text:
{-# LANGUAGE TypeFamilies #-}

import Test.HUnit
import Test.Framework
import Test.Framework.Providers.HUnit

import Numeric.LinearAlgebra.HMatrix
import Numeric.Minimization.QuadProgPP

import Numeric.AD

import Data.SQP

qpTest :: Assertion
qpTest =
  -- From example here:
  -- http://www.mathworks.com/help/optim/ug/quadprog.html#zmw57dd0e69051
  let a = matrix 2 [1, (-1), (-1), 2]
      b = vector [(-2), (-6)]
      e = matrix 2 [ (-1), (-1)
                   , 1, (-2)
                   , (-2), (-1)
                   , 1, 0
                   , 0, 1 ]
      f = vector [2, 2, 3, 0, 0]
      result = solveQuadProg (a, b) Nothing $ Just (e, f)
  in  result @?= Right (vector [0.6666666666666669,
                                1.3333333333333335],
                        -8.22222222222222)

-- sqpTest :: Assertion
-- sqpTest =
--   let a = matrix 2 [1, (-1), (-1), 2]
--       b = vector [(-2), (-6)]
--       e = matrix 2 [ (1), (1)
--                    , (-1), (2)
--                    , (2), (1)
--                    , (-1), 0
--                    , 0, (-1) ]
--       f = vector [(-2), (-2), (-3), 0, 0]
--       affineIneqs x = (e, f)
--       ineqs x = (e #> x) + f
--       cost x = (x `dot` ((0.5 * (a #> x)) + b)) + 0.0
--       problem = Problem
--                 { _cost = cost
--                 , _approxCost = const (a, b, 0)
--                 , _trueIneqs = ineqs
--                 , _approxAffineIneqs = affineIneqs
--                 , _numVariables = 2
--                 , _numIneqs = 5
--                 }
--       (xResult, fResult) = optimize problem (vector [0, 0])
--       (xTrue, fTrue) = (vector [0.6666666666666669,
--                              1.3333333333333335],
--                         -8.22222222222222)
--       xError = norm_2 $ xTrue - xResult
--       fError = abs $ fTrue - fResult
--   in  (xError < 1e-8, fError < 1e-8) @?= (True, True)

getQuadraticApproximation :: ([Double] -> Double)
                          -> ([Double] -> [Double])
                          -> ([Double] -> [[Double]])
                          -> Vector Double
                          -> (Matrix Double, Vector Double, Double)
getQuadraticApproximation fL fL' fL'' x =
  let xL = toList x
      f' = fromList $ fL' xL
      f'' = makeSymmetric $ fromLists $ fL'' xL
      c  = fL xL - (x `dot` f') + 0.5 * (x `dot` (f'' #> x))
  in  ( f'', f' - (f'' #> x), c)

getLinearApproximation :: ([Double] -> [Double])
                       -> ([Double] -> [[Double]])
                       -> Vector Double
                       -> (Matrix Double, Vector Double)
getLinearApproximation fL fL' x =
  let xL = toList x
      f' = fromLists $ fL' xL
  in  (f', vector (fL xL) - (f' #> x))

-- rosenbrockTest :: Assertion
-- rosenbrockTest =
--   -- From: http://www.mathworks.com/help/optim/ug/constrained-nonlinear-optimization-algorithms.html#f26633
--   let f [x, y] = (1 - x)^(2 :: Int) + 100*(y - x^(2 :: Int))^(2 :: Int)
--       cost x = let xL = toList x in f xL
--       convexCost = getQuadraticApproximation f
--       g [x, y] = x^(2 :: Int) + y^(2 :: Int) - 1.5
--       ineqs x = let xL = toList x in vector [g xL]
--       convexIneqs = getLinearApproximation f
--       problem = Problem
--                 { _cost = cost
--                 , _approxCost = convexCost
--                 , _trueIneqs = ineqs
--                 , _approxAffineIneqs = convexIneqs
--                 , _numVariables = 2
--                 , _numIneqs = 1
--                 }
--       (xResult, fResult) = optimize problem (vector [(-1.9), 2.0])
--       (xTrue, fTrue) = (vector [0.9072, 0.8228],
--                         0.00861632761)
--       xError = norm_2 $ xTrue - xResult
--       fError = abs $ fTrue - fResult
--   in  (xError < 1e-4, fError < 1e-4) @?= (True, True)

makeSymmetric :: Matrix Double -> Matrix Double
makeSymmetric m =
  let upperTriIxs = concat
        [ [(r, c) | c <- [(r + 1)..(cols m - 1)]] | r <- [0..(rows m - 1)] ]
      symElems (r, c) = [ ((r, c), m ! r ! c)
                        , ((c, r), m ! r ! c) ]
  in  (diag . takeDiag) m + assoc (size m) 0.0 (concatMap symElems upperTriIxs)

-- From: http://www.mathworks.com/help/optim/ug/nonlinear-equality-and-inequality-constraints.html
nonlinearEqTest :: Assertion
nonlinearEqTest =
  let f [x, y] = exp(x) * (4*x*x + 2*y*y + 4*x*y + 2*y + 1)
      g [x, y] = [negate (x * y) - 10]
      h [x, y] = [(x * x) + y - 1]

      cost = f . toList
      approxCost = getQuadraticApproximation f (grad f) (hessian f)
      ineq = vector . g . toList
      approxIneq = getLinearApproximation g (jacobian g)
      eq = vector . h . toList
      approxEq = getLinearApproximation h (jacobian h)
      problem = Problem
                { _cost = cost
                , _approxCost = approxCost
                , _trueIneqs = ineq
                , _approxAffineIneqs = approxIneq
                , _eqs = eq
                , _approxEqs = approxEq
                , _numVariables = 2
                , _numIneqs = 1
                , _numEqs = 1
                }
      (xResult, fResult) = optimize problem $ vector [(-1.0), 1.0]
      (xTrue, fTrue) = (vector [(-0.7529), 0.4332], 1.5093)
  in  (xResult, fResult) @?= (xTrue, fTrue)

main = defaultMain tests
  where tests = [ -- testCase "qp-test" qpTest
                -- , testCase "sqp-test" sqpTest
                -- , testCase "rosenbrock-test" rosenbrockTest
                  testCase "nonlineareq-test" nonlinearEqTest
                ]

max_stars_repo_name: giogadi/hs-trajopt

max_issues_repo_name: giogadi/hs-trajopt

max_forks_repo_name: giogadi/hs-trajopt


text:
{-# OPTIONS_GHC  -fno-warn-unused-binds -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-missing-signatures #-}
{-# LANGUAGE FlexibleInstances, DeriveFunctor, MultiParamTypeClasses, UndecidableInstances, FlexibleContexts, TypeSynonymInstances #-}

---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------
-- |
-- | Module : First attempt at approximate mean
-- | Creator: Xiao Ling
-- | Created: 11/29/2015
-- | see    : https://github.com/snoyberg/conduit
--            https://gist.github.com/thoughtpolice/3704890
--            http://blog.ezyang.com/2012/01/problem-set-the-codensity-transformation/
--            http://www.janis-voigtlaender.eu/papers/AsymptoticImprovementOfComputationsOverFreeMonads.pdf
--            https://gist.github.com/supki/3776752
-- |
---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------

module Play1129 where

import System.Random

import Data.List
import Data.Random
import Data.Conduit
import qualified Data.Conduit.List as Cl

import Control.Monad.Identity
import Control.Monad.State

import qualified Test.QuickCheck as T

import Core
import Statistics

{-----------------------------------------------------------------------------
   1. make some examples
   2. make
   3. see if cassava can be put into a conduit
   4. output Row type

   see : http://stackoverflow.com/questions/2376981/haskell-types-frustrating-a-simple-average-function
------------------------------------------------------------------------------}

{-----------------------------------------------------------------------------
    II. Approximate Median try ii
    TODO: figure out uniform sampling without knowing length of stream
------------------------------------------------------------------------------}

-- * TODO: make this into an array since we'll be doing (!!) which is O(n) in list
-- * A reservoir is a list, its max allowed size, and its current size

-- * this deign is really bad, cannot gaurantee that currentSize is upToDate with length [a]
-- * dependent types would be good, or we can jsut get rid of current size
-- * or hide it in Core
data Reservoir a = R { run :: [a], size :: Int, currentSize :: Int }
  deriving (Show,Functor,Eq)

instance Num a => T.Arbitrary (Reservoir a) where

  -- * by construction, an arbitrary reservoir is always 1 item away from full
  arbitrary = do
    s <- T.arbitrary :: T.Gen (T.Positive Int)
    let s' = T.getPositive s
    let n = s' - 1 :: Int
    return $ R (replicate n 0) s' n

  shrink (R xs s n) = (\s' -> let n = max 0 (s' - 1) in (R (replicate n 0) s' n)) <$> T.shrink s

-- * constructor an empty resoivor of max size `s`
reservoir :: Int -> Reservoir a
reservoir s = R [] s 0

-- * check if reservoir is full
full :: Reservoir a -> Bool
full (R xs s n) = n >= s

-- * insert item, note the order
(>+) :: Reservoir a -> a -> Reservoir a
(>+) t@(R xs s n) x | full t     = t
                    | otherwise  = R (x:xs) s $ succ n

-- * remove the `i`th item from `r`
-- * if i > size then no item removed
(>-) :: Reservoir a -> Int -> Reservoir a
r@(R xs s n) >- i | i > n     = r
                  | otherwise = undefined      -- * really need to use array here

-- * sample an item `x` from the stream with uniform probability place in reservoir
-- * keep a conter `c` of elements seen
sampl :: Counter -> Sink a (StateT (Reservoir a) RVar) [a]
sampl i = do
  ma <- await
  case ma of
    Nothing -> lift get >>= \r -> return $ run r
    Just a  -> do
      lift $ store i a
      sampl $ succ i

store :: Counter -> a -> StateT (Reservoir a) RVar ()
store i a = do
  r <- get
  case full r of
    False -> put $ r >+ a
    _     -> return ()
      where p = (read . show $ size r :: Float)/i

-- * test m/2 - e*m < rank(y) < m/2 + e*m

-- * quickCheck Resoivor

-- * full resoivor is full
prop_fullReso :: Reservoir Int -> T.Property
prop_fullReso r = (T.==>) (size r > 0) $ full (r>+1) == True

---- * inserting into non-full resoivor resutls in additional element
prop_insertNotFull :: Reservoir Int -> T.Property
prop_insertNotFull r = (T.==>) (size r > 0) $ currentSize (r>+1) == currentSize r + 1

-- * inserting into full resoivor leaves it unchanged
prop_insertFull :: Reservoir Int -> T.Property
prop_insertFull r = (T.==>) (size r > 0) $ currentSize (r'>+1) == currentSize r'
  where r' = r >+ 1

testReservoir :: IO ()
testReservoir = do
  T.quickCheck $ prop_fullReso
  T.quickCheck $ prop_insertFull
  T.quickCheck $ prop_insertNotFull

{-----------------------------------------------------------------------------
    II. Approximate Median Naive
------------------------------------------------------------------------------}

-- * setting problem: can you make this streaming? assuming you know m?
-- * can you do running uniform distribution?

-- * Find approximate eps e-median of list `xs`
-- * Setting: list given before hand for some parameter t
-- * Use : runRVar (naive xs e d) StdRandom
medianNaive :: (Floating a, Ord a) => [a] -> Eps -> Delta -> RVar a
medianNaive xs e d | m <= t    = return $ median xs
             | otherwise = fmap median $ sampl' xs m t
              where
                m = length xs
                t = round $ 7/(e^2) * log (2/d)

-- * use vector here due to (!!)
-- * sample `t` items from the list `[0..m-1]` *with replacement*
-- * Use : runRVar (sample m t) StdRandom
sampl' :: [a] -> Int -> Int -> RVar [a]
sampl' xs m t | m < 0 || t < 0   = return []
          | otherwise        = (\is -> [ xs !! i | i <- is]) <$> mis
            where mis = fmap sort . replicateM t $ uniform 0 (m-1)

{-----------------------------------------------------------------------------
    III. Frequency Moments
------------------------------------------------------------------------------}

{-----------------------------------------------------------------------------
    I. Approximate Counting
------------------------------------------------------------------------------}

-- * Run on stream inputs `xs` for t independent trials for `t = 1/eps`,
-- * and `m` times in parralell, for `m = 1/(e^2 * d)`
-- * and take the median
morris :: Eps -> Delta -> [a] -> RVar Counter
morris e d xs = fmap rmedian . (fmap . fmap) rmean        -- * take median of the means
              $ Cl.sourceList xs $$ count                 -- * stream to counter
              $ replicate m $ replicate t 0               -- * make m counters of length t
              where
                t  = round $ 1/(e^2*d)
                m  = round $ 1/d

-- * Given an m-long list `xs` of lists (each of which is t-lengthed) of counters,
-- * consume the stream and output result
count :: [[Counter]] -> Sink a RVar [[Counter]]
count xxs = (fmap . fmap . fmap) (\x -> 2^(round x) - 1)
            $ Cl.foldM (\xs _ -> incrs xs) xxs

-- * given a list of list of counter `xs` of length `n`,
-- * toss a coin and count
incrs :: [[Counter]] -> RVar [[Counter]]
incrs = sequence . fmap incr where
  incr xs = do
    hs <- sequence $ (\x -> toss . coin $ 0.5^(round x)) <$> xs
    return $ pincr <$> zip hs xs
      where pincr (h,x) = if isHead h then (seq () succ x) else seq () x

-- * Naive solution * --

-- * Run Morris alpha on stream inputs `xs`
morrisA :: [a] -> IO Counter
morrisA xs = flip runRVar StdRandom $ Cl.sourceList xs $$ alpha

-- * Run Morris beta on stream inputs `xs` for `t` independent trials and average
morrisB :: Int -> [a] -> IO Counter
morrisB t =  fmap rmean . replicateM t . morrisA

-- * Naive morris algorithm
-- * Run on stream inputs `xs` for t independent trials for `t = 1/eps`,
-- * and `m` times in parralell, for `m = 1/(e^2 * d)`
-- * and take the median
-- * Problem : `replicateM n` is O(n) time
morris' :: Eps -> Delta -> [a] -> IO Counter
morris' e d = fmap rmedian . replicateM m . morrisB t
  where (t,m) = (round $ 1/(e^2*d), round $ 1/d)

-- * A step in morris Algorithm alpha
alpha :: Sink a RVar Counter
alpha = (\x -> 2^(round x) - 1) <$> Cl.foldM (\x _ -> incr x) 0

-- * Increment a counter `x` with probability 1/2^x
incr :: Counter -> RVar Counter
incr x = do
  h <- toss . coin $ 0.5^(round x)
  return $ if isHead h then (seq () succ x) else seq () x

rmean, rmedian :: (Floating a, Ord a, RealFrac a) => [a] -> Float
rmean   = fromIntegral . round . mean
rmedian = fromIntegral . round . median

---- * Increment a counter `x` with probability 1/2^x
--incr' :: Counter -> RVar Counter
--incr' x = do
--  h <- (\q -> q <= (0.5^(round x) :: Prob)) <$> uniform 0 1
--  return $ if h then (seq () succ x) else seq () x

-- * Test Morris
tMorris :: IO ()
tMorris = undefined

{-----------------------------------------------------------------------------
   X. Finger Exercises

   type Source    m a = ConduitM () a   m ()   -- no meaningful input or return value
   type Conduit a m b = ConduitM a  b   m ()   -- no meaningful return value
   type Sink    a m b = ConduitM a Void m b    -- no meaningful output value
------------------------------------------------------------------------------}

-- * generator
genRs :: Monad m => Source m Int
genRs = Cl.sourceList rs

-- * conduits
cond :: Monad m => Conduit Int m String
cond = Cl.map show

add1 :: Monad m => Conduit Int m Int
add1 = Cl.map (+1)

add2 :: Monad m => Conduit Int m Int
add2 = do
   mx <- await
   case mx of
      Nothing -> return ()
      Just x  -> (yield $ x + 2) >> add2

-- * sinks
sink1 :: Sink Int IO ()
sink1 = Cl.mapM_ $ putStrLn . show

sink2 :: Sink String IO ()
sink2 = Cl.mapM_ putStrLn

-- * accumlate the values in a list
sink3 :: Monad m => [Int] -> Sink Int m [Int]
sink3 xs = do
   mx <- await
   case mx of
      Nothing -> return xs
      Just x  -> sink3 $ xs ++ [x]

-- * accuulate values in a list and increment by 3
add3 :: (Num a, Monad m) => Sink Int m [Int]
add3 = Cl.fold (\xs x -> xs ++ [x+3]) []

-- * some trivial pipes
pp0, pp1, pp2 :: IO ()
pp0 = genRs $$ sink1
pp1 = genRs $$ add1 $= cond $= sink2
pp2 = genRs $$ add2 $= cond $= sink2

-- * map and fold
pp2' :: Identity [Int]
pp2' = genRs $$ add2 $= sink3 []

-- * fold a list
pp3 :: Identity [Int]
pp3 = genRs $$ add3

-- * counts everything
brute :: [a] -> Counter
brute xs = runIdentity $ Cl.sourceList xs $$ Cl.fold (\n _ -> succ n) 0

-- * filter a list
pp4 :: Int -> Identity [Int]
pp4 n = genRs $= Cl.filter (>n) $$ Cl.fold (flip $ (++) . pure) []

-- * map list
pp5 :: Identity [Int]
pp5 = genRs $= Cl.map (+5) $$ sink3 []

-- * test property
tpp2, tpp3, tpp5 :: Bool
tpp2 = runIdentity pp2' == fmap (+2) rs
tpp3 = runIdentity pp3  == fmap (+3) rs
tpp5 = runIdentity pp5  == fmap (+5) rs

-- * test uniformness of list, observe it's relatively uniform
tpp4 :: [Int]
tpp4 = length . runIdentity . pp4 <$> [1..10]

{-----------------------------------------------------------------------------
    X. Test Data
------------------------------------------------------------------------------}

-- * random string of choice
rs :: (Enum a, Num a) => [a]
rs = [1..1000]

{-----------------------------------------------------------------------------
  Depricated
------------------------------------------------------------------------------}

--morrisA' :: Sink a RVar Counter
--morrisA' = Cl.foldM (\x _ -> incr x) 0

--incr :: MonadRandom m => Counter -> m Counter
--incr x = do
--  h <- toss . coin $ 0.5^x
--  return $ if isHead h then (seq () succ x) else seq () x

--morrisB :: Conduit a (StateT Counter RVar) a
--morrisB = do
--  mi <- await
--  case mi of
--    Nothing -> return ()
--    Just i  -> do
--      yield i
--      x <- lift get
--      h <- lift . lift . toss . coin $ 0.5^x
--      let x' = if isHead h then succ x else x
--      lift . put $ seq () x'
--      morrisB

--cap :: Monad m => Sink a m ()
--cap = do
--  mi <- await
--  case mi of
--    Nothing -> return ()
--    _       -> cap

--morrisBs = flip runStateT 0 $ Cl.sourceList [1..10000] $= morrisB $= morrisB $$ cap

-- * Morris Algorithm Beta, repeat step in morris alpha `t` times
--morrisB :: Trials -> Sink a RVar [Counter]
--morrisB t = replicateM t morrisA

--tmorrisB' :: Trials -> [a] -> IO [Counter]
--tmorrisB' t xs = (fmap . fmap) toN . flip runRVar StdRandom $ Cl.sourceList xs $$ morrisB t


text:
module HankelTransform where

import           Control.Parallel.Strategies
import           Data.Array.Repa                        as R
import           Data.Complex                           as C
import           Data.List                              as L
import           Data.Vector.Storable                   as VS
import           Data.Vector.Unboxed                    as VU
import           Graphics.Rendering.Chart.Backend.Cairo
import           Graphics.Rendering.Chart.Easy
import           Numeric.GSL.Special.Bessel
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.BLAS
import           Utils.List
import           Utils.Parallel

main = do
  args@(nStr:deltaXStr:deltaRStr:radiusXStr:radiusRStr:_) <- getArgs
  let n = read nStr :: Int
      deltaX = read deltaXStr :: Double
      deltaR = read deltaRStr :: Double
      radiusX = read radiusXStr :: Double
      radiusR = read radiusRStr :: Double
      folderPath = "output/test/HankelTransform"
  createDirectoryIfMissing True folderPath
  let numR = round $ radiusR / deltaR :: Int
      numX = round $ radiusX / deltaX :: Int
      func r =
        if r >= 1
          then 1 / r
          else 1 / r
      fr =
        parMap
          rdeepseq
          (\r -> func r * r)
          [fromIntegral r * deltaR | r <- [1 .. numR]]
      besselArrayR =
        fromFunction (Z :. numX :. numR) $ \(Z :. x :. r) ->
          deltaR *
          bessel_Jn
            (fromIntegral n)
            ((fromIntegral x + 1) * deltaX * (fromIntegral r + 1) * deltaR)
  besselArray <- computeUnboxedP besselArrayR
  let besselVec = VU.convert . toUnboxed $ besselArray
      frVec = VS.fromList fr
  xs <- VS.toList <$> gemmBLAS numX 1 numR besselVec frVec
  toFile def (folderPath </> "HankelTransform.png") $ do
    layout_title .=
      printf "N = %d RadiusR = %.2f RadiusX = %.2f" n radiusR radiusX
    plot (line "" [L.zip [fromIntegral i * deltaX | i <- [1 .. numX]] xs])

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE TypeSynonymInstances #-}
module Main where

import           Control.Applicative     ((<$>))
import           Data.Complex
import           Data.Vector
import           Prelude                 hiding (length, map, maximum, sum,
                                          zipWith, (++))
import           Test.QuickCheck.Monadic
import           Test.Tasty
import           Test.Tasty.QuickCheck   hiding (generate)

import           Numeric.FFT

main :: IO ()
main = defaultMain tests

tests :: TestTree
tests = testGroup "Tests"
        [ testProperty "FFT vs. DFT" prop_dft_vs_fft
        , testProperty "FFT/IFFT round-trip" prop_ifft
        ]

-- Clean up number display.
defuzz :: Vector (Complex Double) -> Vector (Complex Double)
defuzz = map (\(r :+ i) -> df r :+ df i)
  where df x = if abs x < 1.0E-6 then 0 else x

-- Check FFT against DFT.
check :: Vector (Complex Double) -> IO (Double, Vector (Complex Double))
check v = do
  tst <- fft v
  let diff = defuzz $ zipWith (-) tst (basicDFT v)
  return (maximum $ map magnitude diff, diff)

-- QuickCheck property for FFT vs. DFT testing.
prop_dft_vs_fft :: Property
prop_dft_vs_fft = monadicIO $ do
  v <- pick arbitrary
  chk <- run $ check v
  assert $ fst chk < 1.0E-6

-- QuickCheck property for inverse FFT round-trip testing.
prop_ifft :: Property
prop_ifft = monadicIO $ do
  v <- pick arbitrary
  fwd <- run $ fft v
  bwd <- run $ ifft fwd
  let diff = zipWith (-) v bwd
  assert $ maximum (map magnitude diff) < 1.0E-6

-- Non-zero length arbitrary vectors.
instance Arbitrary (Vector (Complex Double)) where
  arbitrary = fromList <$> listOf1 arbitrary

-- Roots of unity.
omega :: Int -> Complex Double
omega n = cis (2 * pi / fromIntegral n)

-- Naïve DFT.
basicDFT :: Vector (Complex Double) -> Vector (Complex Double)
basicDFT h = generate n doone
  where n = length h
        w = omega n
        doone i = sum $ zipWith (*) h $ generate n (\k -> w^^(i*k))

max_stars_repo_name: ian-ross/arb-fft

max_issues_repo_name: ian-ross/arb-fft

max_forks_repo_name: ian-ross/arb-fft


text:
module Perceptron where

import Numeric.LinearAlgebra

v00 = vector [0.0,0.0]
v01 = vector [0.0,1.0]
v10 = vector [1.0,0.0]
v11 = vector [1.0,1.0]

p_and :: Vector Double -> Double
p_and x
  | tmp <= 0.0 = 0.0
  | otherwise = 1.0
  where
    tmp = w <.> x + b
    b = -0.7
    w = vector [0.5,0.5]

p_nand :: Vector Double -> Double
p_nand x
  | tmp <= 0.0 = 0.0
  | otherwise = 1.0
  where
    tmp = w <.> x + b
    b = 0.7
    w = vector [-0.5,-0.5]

p_or :: Vector Double -> Double
p_or x
  | tmp <= 0.0 = 0.0
  | otherwise = 1.0
  where
    tmp = w <.> x + b
    b = -0.2
    w = vector [0.5,0.5]

p_xor :: Vector Double -> Double
p_xor x = p_and $ vector [s1,s2]
  where
    s1 = p_nand x
    s2 = p_or x

max_stars_repo_name: chupaaaaaaan/nn-with-haskell

max_issues_repo_name: chupaaaaaaan/nn-with-haskell

max_forks_repo_name: chupaaaaaaan/nn-with-haskell


text:
{-# LANGUAGE CPP, TypeSynonymInstances, FlexibleInstances #-}

module MyVectorType (module MyVectorType
#ifdef VECTOR_HMATRIX
 , Matrix.toLists, Matrix.fromLists
 , Matrix.trans
 , Vector.toList, Vector.fromList
 , Vector.dim
 , Numeric.Container.vXm, Numeric.Container.sumElements
#endif
)
where

#ifdef VECTOR_HMATRIX
import qualified Data.Packed.Vector as Vector
import qualified Data.Packed.Matrix as Matrix
import qualified Numeric.Container
import Numeric.Container (Element, Container)
import Numeric.LinearAlgebra () -- instances
type Vector = Vector.Vector
type Matrix = Matrix.Matrix
instance Ord (Matrix Double) where
    compare a b = compare (Matrix.toLists a) (Matrix.toLists b)
cmap :: (Element b, Container c e) => (e -> b) -> c e -> c b
cmap a  = Numeric.Container.cmap a

add :: Container c e => c e -> c e -> c e
add a = Numeric.Container.add a

transFix :: Matrix Double -> Matrix Double
transFix = id

-- concat :: (Storable a) => [Vector a] -> Vector a
concat :: [Vector Double] -> Vector Double
concat = Vector.join
#endif

#ifdef VECTOR_VECTOR
import Data.Vector.Binary () -- instances
import Data.List (transpose)
#ifdef VECTOR_DDOT
import Foreign.C.Types
import Foreign
import qualified Data.Vector.Storable as V
import qualified Data.Vector.Storable.Internal as VI
#else
import qualified Data.Vector as V
#endif

type Vector a = V.Vector a
type Matrix a = [Vector a]

#ifdef VECTOR_DDOT
foreign import ccall "cblas_ddot" cblas_ddot :: CInt -> Ptr CDouble -> CInt -> Ptr CDouble -> CInt -> CDouble
vDOTv :: Vector Double -> Vector Double -> Double
vDOTv v1 v2 = let (fptr'1,len1) = V.unsafeToForeignPtr0 (V.unsafeCast v1)
                  (fptr'2,len2) = V.unsafeToForeignPtr0 (V.unsafeCast v2)
                  ptr'v1 = VI.getPtr fptr'1
                  ptr'v2 = VI.getPtr fptr'2
                  r = (cblas_ddot (fromIntegral len1) ptr'v1 1 ptr'v2 1)
              in realToFrac r
vXm :: Vector Double -> Matrix Double -> Vector Double
vXm v m = V.fromList $ map (vDOTv v) m
#else
vXm :: Vector Double -> Matrix Double -> Vector Double
vXm v m = V.fromList $ map (\ vm -> V.sum $ V.zipWith (*) v vm) m
#endif

add :: Vector Double -> Vector Double -> Vector Double
add = V.zipWith (+)

cmap :: (Double -> Double) -> Vector Double -> Vector Double
cmap = V.map

toLists :: Matrix Double -> [[Double]]
toLists = map V.toList
fromLists :: [[Double]] -> Matrix Double
fromLists = map V.fromList
trans :: Matrix Double -> Matrix Double
trans = fromLists . transpose . toLists

transFix :: Matrix Double -> Matrix Double
transFix = trans

toList :: Vector Double -> [Double]
toList = V.toList
fromList :: [Double] -> Vector Double
fromList = V.fromList

sumElements :: Vector Double -> Double
sumElements = V.sum

dim :: Vector Double -> Int
dim = V.length

concat :: [Vector a] -> Vector a
concat = V.concat
#endif

-- #ifdef VECTOR_DPH
-- import qualified Prelude
-- import Prelude ((.), ($))
-- -- import Data.Array.Parallel
-- -- import Data.Array.Parallel.Prelude
-- import Data.Array.Parallel.Prelude.Double
-- import Data.Array.Parallel.Prelude.Int (Int)
-- import Data.Array.Parallel.Unlifted as V
-- import Data.List (transpose)
--
-- type Vector a = V.Array a
-- type Matrix a = [Vector a]
--
-- vXm :: Vector Double -> Matrix Double -> Vector Double
-- vXm v m = Prelude.map (\ vm -> V.sum $ V.zipWith (*) v vm) m
--
-- add :: Vector Double -> Vector Double -> Vector Double
-- add = V.zipWith (+)
--
-- cmap :: (Double -> Double) -> Vector Double -> Vector Double
-- cmap = V.map
--
-- toLists :: Matrix Double -> [[Double]]
-- toLists = V.toList . V.map V.toList
-- fromLists :: [[Double]] -> Matrix Double
-- fromLists = V.map V.fromList . V.fromList
-- trans :: Matrix Double -> Matrix Double
-- trans = fromLists . transpose . toLists
--
-- transFix :: Matrix Double -> Matrix Double
-- transFix = trans
--
-- toList :: Vector Double -> [Double]
-- toList = V.toList
-- fromList :: [Double] -> Vector Double
-- fromList = V.fromList
--
-- sumElements :: Vector Double -> Double
-- sumElements = V.sum
--
-- dim :: Vector Double -> Int
-- dim = V.length
--
-- #endif

max_stars_repo_name: Tener/deeplearning-thesis

max_issues_repo_name: Tener/deeplearning-thesis

max_forks_repo_name: Tener/deeplearning-thesis


text:
{-# LANGUAGE CPP #-}
-- |
-- Module: Optics.At.Core
-- Description: Optics for 'Map' and 'Set'-like containers.
--
-- This module provides optics for 'Map' and 'Set'-like containers, including an
-- 'AffineTraversal' to traverse a key in a map or an element of a sequence:
--
-- >>> preview (ix 1) ['a','b','c']
-- Just 'b'
--
-- a 'Lens' to get, set or delete a key in a map:
--
-- >>> set (at 0) (Just 'b') (Map.fromList [(0, 'a')])
-- fromList [(0,'b')]
--
-- and a 'Lens' to insert or remove an element of a set:
--
-- >>> IntSet.fromList [1,2,3,4] & contains 3 .~ False
-- fromList [1,2,4]
--
-- The @Optics.At@ module from @optics-extra@ provides additional instances of
-- the classes defined here.
--
module Optics.At.Core
  (
    -- * Type families
    Index
  , IxValue

    -- * Ixed
  , Ixed(..)
  , ixAt

    -- * At
  , At(..)
  , at'
  , sans

  -- * Contains
  , Contains(..)
  ) where

import Data.Array.IArray as Array
import Data.Array.Unboxed
import Data.Complex
import Data.Functor.Identity
import Data.IntMap as IntMap
import Data.IntSet as IntSet
import Data.List.NonEmpty as NonEmpty
import Data.Map as Map
import Data.Sequence as Seq
import Data.Set as Set
import Data.Tree

import Data.Maybe.Optics
import Optics.AffineTraversal
import Optics.Iso
import Optics.Lens
import Optics.Optic
import Optics.Setter

-- | Type family that takes a key-value container type and returns the type of
-- keys (indices) into the container, for example @'Index' ('Map' k a) ~ k@.
-- This is shared by 'Ixed', 'At' and 'Contains'.
type family Index (s :: *) :: *
type instance Index (e -> a) = e
type instance Index IntSet = Int
type instance Index (Set a) = a
type instance Index [a] = Int
type instance Index (NonEmpty a) = Int
type instance Index (Seq a) = Int
type instance Index (a,b) = Int
type instance Index (a,b,c) = Int
type instance Index (a,b,c,d) = Int
type instance Index (a,b,c,d,e) = Int
type instance Index (a,b,c,d,e,f) = Int
type instance Index (a,b,c,d,e,f,g) = Int
type instance Index (a,b,c,d,e,f,g,h) = Int
type instance Index (a,b,c,d,e,f,g,h,i) = Int
type instance Index (IntMap a) = Int
type instance Index (Map k a) = k
type instance Index (Array.Array i e) = i
type instance Index (UArray i e) = i
type instance Index (Complex a) = Int
type instance Index (Identity a) = ()
type instance Index (Maybe a) = ()
type instance Index (Tree a) = [Int]

-- | This class provides a simple 'Lens' that lets you view (and modify)
-- information about whether or not a container contains a given 'Index'.
-- Instances are provided for 'Set'-like containers only.
class Contains m where
  -- |
  -- >>> IntSet.fromList [1,2,3,4] ^. contains 3
  -- True
  --
  -- >>> IntSet.fromList [1,2,3,4] ^. contains 5
  -- False
  --
  -- >>> IntSet.fromList [1,2,3,4] & contains 3 .~ False
  -- fromList [1,2,4]
  contains :: Index m -> Lens' m Bool

instance Contains IntSet where
  contains k = lensVL $ \f s -> f (IntSet.member k s) <&> \b ->
    if b then IntSet.insert k s else IntSet.delete k s
  {-# INLINE contains #-}

instance Ord a => Contains (Set a) where
  contains k = lensVL $ \f s -> f (Set.member k s) <&> \b ->
    if b then Set.insert k s else Set.delete k s
  {-# INLINE contains #-}

-- | Type family that takes a key-value container type and returns the type of
-- values stored in the container, for example @'IxValue' ('Map' k a) ~ a@. This
-- is shared by both 'Ixed' and 'At'.
type family IxValue (m :: *) :: *

-- | Provides a simple 'AffineTraversal' lets you traverse the value at a given
-- key in a 'Map' or element at an ordinal position in a list or 'Seq'.
class Ixed m where
  -- | /NB:/ Setting the value of this 'AffineTraversal' will only set the value
  -- in 'at' if it is already present.
  --
  -- If you want to be able to insert /missing/ values, you want 'at'.
  --
  -- >>> [1,2,3,4] & ix 2 %~ (*10)
  -- [1,2,30,4]
  --
  -- >>> "abcd" & ix 2 .~ 'e'
  -- "abed"
  --
  -- >>> "abcd" ^? ix 2
  -- Just 'c'
  --
  -- >>> [] ^? ix 2
  -- Nothing
  ix :: Index m -> AffineTraversal' m (IxValue m)
  default ix :: At m => Index m -> AffineTraversal' m (IxValue m)
  ix = ixAt
  {-# INLINE ix #-}

-- | A definition of 'ix' for types with an 'At' instance. This is the default
-- if you don't specify a definition for 'ix'.
ixAt :: At m => Index m -> AffineTraversal' m (IxValue m)
ixAt = \i -> at i % _Just
{-# INLINE ixAt #-}

type instance IxValue (e -> a) = a
instance Eq e => Ixed (e -> a) where
  ix e = atraversalVL $ \_ p f -> p (f e) <&> \a e' -> if e == e' then a else f e'
  {-# INLINE ix #-}

type instance IxValue (Maybe a) = a
instance Ixed (Maybe a) where
  ix () = castOptic @An_AffineTraversal _Just
  {-# INLINE ix #-}

type instance IxValue [a] = a
instance Ixed [a] where
  ix k = atraversalVL (ixListVL k)
  {-# INLINE ix #-}

type instance IxValue (NonEmpty a) = a
instance Ixed (NonEmpty a) where
  ix k = atraversalVL $ \point f xs0 ->
    if k < 0
    then point xs0
    else let go (a:|as) 0 = f a <&> (:|as)
             go (a:|as) i = (a:|) <$> ixListVL (i - 1) point f as
         in go xs0 k
  {-# INLINE ix #-}

type instance IxValue (Identity a) = a
instance Ixed (Identity a) where
  ix () = atraversalVL $ \_ f (Identity a) -> Identity <$> f a
  {-# INLINE ix #-}

type instance IxValue (Tree a) = a
instance Ixed (Tree a) where
  ix xs0 = atraversalVL $ \point f ->
    let go [] (Node a as) = f a <&> \a' -> Node a' as
        go (i:is) t@(Node a as)
          | i < 0     = point t
          | otherwise = Node a <$> ixListVL i point (go is) as
    in go xs0
  {-# INLINE ix #-}

type instance IxValue (Seq a) = a
instance Ixed (Seq a) where
  ix i = atraversalVL $ \point f m ->
    if 0 <= i && i < Seq.length m
    then f (Seq.index m i) <&> \a -> Seq.update i a m
    else point m
  {-# INLINE ix #-}

type instance IxValue (IntMap a) = a
-- Default implementation uses IntMap.alterF
instance Ixed (IntMap a)

type instance IxValue (Map k a) = a
-- Default implementation uses Map.alterF
instance Ord k => Ixed (Map k a)

type instance IxValue (Set k) = ()
instance Ord k => Ixed (Set k) where
  ix k = atraversalVL $ \point f m ->
    if Set.member k m
    then f () <&> \() -> Set.insert k m
    else point m
  {-# INLINE ix #-}

type instance IxValue IntSet = ()
instance Ixed IntSet where
  ix k = atraversalVL $ \point f m ->
    if IntSet.member k m
    then f () <&> \() -> IntSet.insert k m
    else point m
  {-# INLINE ix #-}

type instance IxValue (Array.Array i e) = e
-- |
-- @
-- arr 'Array.!' i ≡ arr 'Optics.Operators.^.' 'ix' i
-- arr '//' [(i,e)] ≡ 'ix' i 'Optics.Operators..~' e '$' arr
-- @
instance Ix i => Ixed (Array.Array i e) where
  ix i = atraversalVL $ \point f arr ->
    if inRange (bounds arr) i
    then f (arr Array.! i) <&> \e -> arr Array.// [(i,e)]
    else point arr
  {-# INLINE ix #-}

type instance IxValue (UArray i e) = e
-- |
-- @
-- arr 'Array.!' i ≡ arr 'Optics.Operators.^.' 'ix' i
-- arr '//' [(i,e)] ≡ 'ix' i 'Optics.Operators..~' e '$' arr
-- @
instance (IArray UArray e, Ix i) => Ixed (UArray i e) where
  ix i = atraversalVL $ \point f arr ->
    if inRange (bounds arr) i
    then f (arr Array.! i) <&> \e -> arr Array.// [(i,e)]
    else point arr
  {-# INLINE ix #-}

-- | @'ix' :: 'Int' -> 'AffineTraversal'' (a, a) a@
type instance IxValue (a0, a2) = a0
instance (a0 ~ a1) => Ixed (a0, a1) where
  ix i = atraversalVL $ \point f ~s@(a0, a1) ->
    case i of
      0 -> (,a1) <$> f a0
      1 -> (a0,) <$> f a1
      _ -> point s

-- | @'ix' :: 'Int' -> 'AffineTraversal'' (a, a, a) a@
type instance IxValue (a0, a1, a2) = a0
instance (a0 ~ a1, a0 ~ a2) => Ixed (a0, a1, a2) where
  ix i = atraversalVL $ \point f ~s@(a0, a1, a2) ->
    case i of
      0 -> (,a1,a2) <$> f a0
      1 -> (a0,,a2) <$> f a1
      2 -> (a0,a1,) <$> f a2
      _ -> point s

-- | @'ix' :: 'Int' -> 'AffineTraversal'' (a, a, a, a) a@
type instance IxValue (a0, a1, a2, a3) = a0
instance (a0 ~ a1, a0 ~ a2, a0 ~ a3) => Ixed (a0, a1, a2, a3) where
  ix i = atraversalVL $ \point f ~s@(a0, a1, a2, a3) ->
    case i of
      0 -> (,a1,a2,a3) <$> f a0
      1 -> (a0,,a2,a3) <$> f a1
      2 -> (a0,a1,,a3) <$> f a2
      3 -> (a0,a1,a2,) <$> f a3
      _ -> point s

-- | @'ix' :: 'Int' -> 'AffineTraversal'' (a, a, a, a, a) a@
type instance IxValue (a0, a1, a2, a3, a4) = a0
instance (a0 ~ a1, a0 ~ a2, a0 ~ a3, a0 ~ a4) => Ixed (a0, a1, a2, a3, a4) where
  ix i = atraversalVL $ \point f ~s@(a0, a1, a2, a3, a4) ->
    case i of
      0 -> (,a1,a2,a3,a4) <$> f a0
      1 -> (a0,,a2,a3,a4) <$> f a1
      2 -> (a0,a1,,a3,a4) <$> f a2
      3 -> (a0,a1,a2,,a4) <$> f a3
      4 -> (a0,a1,a2,a3,) <$> f a4
      _ -> point s

-- | @'ix' :: 'Int' -> 'AffineTraversal'' (a, a, a, a, a, a) a@
type instance IxValue (a0, a1, a2, a3, a4, a5) = a0
instance
  (a0 ~ a1, a0 ~ a2, a0 ~ a3, a0 ~ a4, a0 ~ a5
  ) => Ixed (a0, a1, a2, a3, a4, a5) where
  ix i = atraversalVL $ \point f ~s@(a0, a1, a2, a3, a4, a5) ->
    case i of
      0 -> (,a1,a2,a3,a4,a5) <$> f a0
      1 -> (a0,,a2,a3,a4,a5) <$> f a1
      2 -> (a0,a1,,a3,a4,a5) <$> f a2
      3 -> (a0,a1,a2,,a4,a5) <$> f a3
      4 -> (a0,a1,a2,a3,,a5) <$> f a4
      5 -> (a0,a1,a2,a3,a4,) <$> f a5
      _ -> point s

-- | @'ix' :: 'Int' -> 'AffineTraversal'' (a, a, a, a, a, a, a) a@
type instance IxValue (a0, a1, a2, a3, a4, a5, a6) = a0
instance
  (a0 ~ a1, a0 ~ a2, a0 ~ a3, a0 ~ a4, a0 ~ a5, a0 ~ a6
  ) => Ixed (a0, a1, a2, a3, a4, a5, a6) where
  ix i = atraversalVL $ \point f ~s@(a0, a1, a2, a3, a4, a5, a6) ->
    case i of
      0 -> (,a1,a2,a3,a4,a5,a6) <$> f a0
      1 -> (a0,,a2,a3,a4,a5,a6) <$> f a1
      2 -> (a0,a1,,a3,a4,a5,a6) <$> f a2
      3 -> (a0,a1,a2,,a4,a5,a6) <$> f a3
      4 -> (a0,a1,a2,a3,,a5,a6) <$> f a4
      5 -> (a0,a1,a2,a3,a4,,a6) <$> f a5
      6 -> (a0,a1,a2,a3,a4,a5,) <$> f a6
      _ -> point s

-- | @'ix' :: 'Int' -> 'AffineTraversal'' (a, a, a, a, a, a, a, a) a@
type instance IxValue (a0, a1, a2, a3, a4, a5, a6, a7) = a0
instance
  (a0 ~ a1, a0 ~ a2, a0 ~ a3, a0 ~ a4, a0 ~ a5, a0 ~ a6, a0 ~ a7
  ) => Ixed (a0, a1, a2, a3, a4, a5, a6, a7) where
  ix i = atraversalVL $ \point f ~s@(a0, a1, a2, a3, a4, a5, a6, a7) ->
    case i of
      0 -> (,a1,a2,a3,a4,a5,a6,a7) <$> f a0
      1 -> (a0,,a2,a3,a4,a5,a6,a7) <$> f a1
      2 -> (a0,a1,,a3,a4,a5,a6,a7) <$> f a2
      3 -> (a0,a1,a2,,a4,a5,a6,a7) <$> f a3
      4 -> (a0,a1,a2,a3,,a5,a6,a7) <$> f a4
      5 -> (a0,a1,a2,a3,a4,,a6,a7) <$> f a5
      6 -> (a0,a1,a2,a3,a4,a5,,a7) <$> f a6
      7 -> (a0,a1,a2,a3,a4,a5,a6,) <$> f a7
      _ -> point s

-- | @'ix' :: 'Int' -> 'AffineTraversal'' (a, a, a, a, a, a, a, a, a) a@
type instance IxValue (a0, a1, a2, a3, a4, a5, a6, a7, a8) = a0
instance
  (a0 ~ a1, a0 ~ a2, a0 ~ a3, a0 ~ a4, a0 ~ a5, a0 ~ a6, a0 ~ a7, a0 ~ a8
  ) => Ixed (a0, a1, a2, a3, a4, a5, a6, a7, a8) where
  ix i = atraversalVL $ \point f ~s@(a0, a1, a2, a3, a4, a5, a6, a7, a8) ->
    case i of
      0 -> (,a1,a2,a3,a4,a5,a6,a7,a8) <$> f a0
      1 -> (a0,,a2,a3,a4,a5,a6,a7,a8) <$> f a1
      2 -> (a0,a1,,a3,a4,a5,a6,a7,a8) <$> f a2
      3 -> (a0,a1,a2,,a4,a5,a6,a7,a8) <$> f a3
      4 -> (a0,a1,a2,a3,,a5,a6,a7,a8) <$> f a4
      5 -> (a0,a1,a2,a3,a4,,a6,a7,a8) <$> f a5
      6 -> (a0,a1,a2,a3,a4,a5,,a7,a8) <$> f a6
      7 -> (a0,a1,a2,a3,a4,a5,a6,,a8) <$> f a7
      8 -> (a0,a1,a2,a3,a4,a5,a6,a7,) <$> f a8
      _ -> point s

-- | 'At' provides a 'Lens' that can be used to read, write or delete the value
-- associated with a key in a 'Map'-like container on an ad hoc basis.
--
-- An instance of 'At' should satisfy:
--
-- @
-- 'ix' k ≡ 'at' k '%' '_Just'
-- @
class Ixed m => At m where
  -- |
  -- >>> Map.fromList [(1,"world")] ^. at 1
  -- Just "world"
  --
  -- >>> at 1 ?~ "hello" $ Map.empty
  -- fromList [(1,"hello")]
  --
  -- /Note:/ Usage of this function might introduce space leaks if you're not
  -- careful to make sure that values put inside the 'Just' constructor are
  -- evaluated. To force the values and avoid such leaks, use 'at'' instead.
  --
  -- /Note:/ 'Map'-like containers form a reasonable instance, but not
  -- 'Array'-like ones, where you cannot satisfy the 'Lens' laws.
  at :: Index m -> Lens' m (Maybe (IxValue m))

-- | Version of 'at' strict in the value inside the `Just` constructor.
--
-- Example:
--
-- >>> (at () .~ Just (error "oops") $ Nothing) `seq` ()
-- ()
--
-- >>> (at' () .~ Just (error "oops") $ Nothing) `seq` ()
-- *** Exception: oops
-- ...
--
-- >>> view (at ()) (Just $ error "oops") `seq` ()
-- ()
--
-- >>> view (at' ()) (Just $ error "oops") `seq` ()
-- *** Exception: oops
-- ...
--
-- It also works as expected for other data structures:
--
-- >>> (at 1 .~ Just (error "oops") $ Map.empty) `seq` ()
-- ()
--
-- >>> (at' 1 .~ Just (error "oops") $ Map.empty) `seq` ()
-- *** Exception: oops
-- ...
at' :: At m => Index m -> Lens' m (Maybe (IxValue m))
at' k = at k % iso f f
  where
    f = \case
      Just !x -> Just x
      Nothing -> Nothing
{-# INLINE at' #-}

-- | Delete the value associated with a key in a 'Map'-like container
--
-- @
-- 'sans' k = 'at' k 'Optics.Operators..~' Nothing
-- @
sans :: At m => Index m -> m -> m
sans k = set (at k) Nothing
{-# INLINE sans #-}

instance At (Maybe a) where
  at () = lensVL id
  {-# INLINE at #-}

instance At (IntMap a) where
#if MIN_VERSION_containers(0,5,8)
  at k = lensVL $ \f -> IntMap.alterF f k
#else
  at k = lensVL $ \f m ->
    let mv = IntMap.lookup k m
    in f mv <&> \r -> case r of
      Nothing -> maybe m (const (IntMap.delete k m)) mv
      Just v' -> IntMap.insert k v' m
#endif
  {-# INLINE at #-}

instance Ord k => At (Map k a) where
#if MIN_VERSION_containers(0,5,8)
  at k = lensVL $ \f -> Map.alterF f k
#else
  at k = lensVL $ \f m ->
    let mv = Map.lookup k m
    in f mv <&> \r -> case r of
      Nothing -> maybe m (const (Map.delete k m)) mv
      Just v' -> Map.insert k v' m
#endif
  {-# INLINE at #-}

instance At IntSet where
  at k = lensVL $ \f m ->
    let mv = if IntSet.member k m
             then Just ()
             else Nothing
    in f mv <&> \r -> case r of
      Nothing -> maybe m (const (IntSet.delete k m)) mv
      Just () -> IntSet.insert k m
  {-# INLINE at #-}

instance Ord k => At (Set k) where
  at k = lensVL $ \f m ->
    let mv = if Set.member k m
             then Just ()
             else Nothing
    in f mv <&> \r -> case r of
      Nothing -> maybe m (const (Set.delete k m)) mv
      Just () -> Set.insert k m
  {-# INLINE at #-}

----------------------------------------
-- Internal

ixListVL :: Int -> AffineTraversalVL' [a] a
ixListVL k point f xs0 =
  if k < 0
  then point xs0
  else let go [] _ = point []
           go (a:as) 0 = f a <&> (:as)
           go (a:as) i = (a:) <$> (go as $! i - 1)
       in go xs0 k
{-# INLINE ixListVL #-}

-- $setup
-- >>> import Optics.Core

max_stars_repo_name: LiangJianyi/HaskellLensOptics

max_issues_repo_name: LiangJianyi/HaskellLensOptics

max_forks_repo_name: LiangJianyi/HaskellLensOptics


text:
{-# LANGUAGE FlexibleContexts #-}
module Network (initNet, feedForward, sgd, forwardPass, backprop, updateMiniBatch) where

import Control.Monad
import Data.List.Split
import Numeric.LinearAlgebra as L

import MnistLoader
import Aux

feedForward :: Network -> Matrix Float -> Matrix Float
feedForward n a = foldl f a $ zip (weights n) (biases n)
  where
    f a' (w, b) = sigmoid ((w L.<> a') + b)

sgd :: Network
  -> [(Matrix Float, Matrix Float)] -- traning data
  -> [(Matrix Float, Matrix Float)] -- Evaluation data
  -> Int -- Number of epochs
  -> Int -- Current Epoc
  -> Int -- Mini batch size
  -> Float -- Eta
  -> IO Network
sgd n exemplars testData maxEpoch epoch batchSize eta = do
  exemplars' <- shuffle exemplars

  let batches = chunksOf batchSize exemplars'
  let newNet = foldl (\accNet batch -> updateMiniBatch accNet batch eta) n batches
  let correct = evaluate newNet testData
  putStrLn $ "Epoch " ++ (show (epoch + 1)) ++ " Eval: " ++ (show correct) ++ " / " ++ (show $ length testData)

  if epoch == maxEpoch
    then return newNet
    else sgd newNet exemplars' testData maxEpoch (epoch + 1) batchSize eta

updateMiniBatch :: Network -> [(Matrix Float, Matrix Float)] -> Float -> Network
updateMiniBatch n l@((x, y) : rest) eta = let
    initAcc = backprop n x y
    nablaSums = foldl f initAcc rest
    nabla_b = map fst nablaSums
    nabla_w = map snd nablaSums
  in
    Network {
        sizes = sizes n
      , weights = zipWith wavg (weights n) nabla_w
      , biases = zipWith wavg (biases n) nabla_b
    }
  where
    f nablaSum (x, y) = let
        bp = backprop n x y
      in
        zipWith add nablaSum bp
    add (nb, nw) (nb', nw') = (nb+nb', nw+nw')
    wavg :: Matrix Float -> Matrix Float -> Matrix Float
    wavg m1 m2 = let
        nt = (1><1) [fromIntegral $ length l]
        eta' = (1><1) [eta]
      in m1 - (eta' / nt) * m2

-- Network -> x -> y -> (nabla_b, nabla_w)
backprop :: Network -> Matrix Float -> Matrix Float -> [(Matrix Float, Matrix Float)]
backprop n x y = let
    -- Reversed list of forward passes
    ((a:a':as), (z:zs)) = forwardPass n x
    revW = foldl (\acc w -> w:acc) [] (weights n)
    fp = zip as zs
    zippedFp = zip fp revW

    -- Handle last layer manually
    delta = (cost_derivative a y) * (sigmoid' z)
    nabla_b = delta
    nabla_w = delta L.<> (tr' a')
    res = foldl f [(nabla_b, nabla_w)] zippedFp
  in
    res
  where
    f (acc@((delta, _):_)) ((a, z), w) = let
        sp = sigmoid' z
        delta' = ((tr' w) L.<> delta) * sp
        nabla_b = delta'
        nabla_w = delta' L.<> (tr' a)
      in
        (nabla_b, nabla_w) : acc

-- Network -> activation -> (activations, z vectors)
forwardPass :: Network -> Matrix Float -> ([Matrix Float], [Matrix Float])
forwardPass n a = let
  pairs = zip (weights n) (biases n)
  fp = foldl f [(a, a)] pairs
  in
    (map fst fp, init $ map snd fp)
    where
      f (l@((a, _):_)) (w, b) = let
          z = (w L.<> a) + b
          activation = sigmoid z
        in (activation, z):l

evaluate :: Network -> [(Matrix Float, Matrix Float)] -> Int
evaluate net exemplars = sum $ map f exemplars
  where
    f (x, y) = let
        approx = argMax $ feedForward net x
        actual = argMax y
      in
        if approx == actual
          then 1
          else 0

cost_derivative output_activations y = output_activations - y

-- Miscellaneous
sigmoid z = 1.0 / (1.0 + (exp $ -z))

sigmoid' :: Matrix Float -> Matrix Float
sigmoid' z = (sigmoid z) * (1 - (sigmoid z))

argMax :: Matrix Float -> Int
argMax m = L.maxIndex $ L.flatten m

-- assuming row-major order
data Network = Network { sizes :: [Int]
                       , weights :: [Matrix Float]
                       , biases :: [Matrix Float]
                       } deriving (Show)

initNet :: [Int] -> IO (Network)
initNet sizes = do
  ws <- genWeights sizes
  bs <- genBiases sizes
  return Network { sizes = sizes
                  , weights = ws
                  , biases = bs
                  }
  where
    genWeights [s, s'] =
      do
        m <- L.randn s' s
        return [L.cmap realToFrac m]
    genWeights (s:s':ss) =
      do
        rest <- genWeights (s' : ss)
        m <- L.randn s' s
        return (L.cmap realToFrac m : rest)
    genWeights _ = error "Not enough layers"
    genBiases [s, s'] =
      do
        m <- L.randn s' 1
        return $ [L.cmap realToFrac m]
    genBiases (s:s':ss) =
      do
        rest <- genBiases (s' : ss)
        m <- L.randn s' 1
        return (L.cmap realToFrac m : rest)
    genBiases _ = error "Not enough layers"

max_stars_repo_name: madsbuch/haskell-nn

max_issues_repo_name: madsbuch/haskell-nn

max_forks_repo_name: madsbuch/haskell-nn


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

module Roc (testRoc) where

import Data.Csv ((.:), FromNamedRecord(..), decodeByName)
import Codec.Compression.Lzma (decompress)
import Test.Tasty
import Test.Tasty.HUnit
import Test.Tasty.QuickCheck
import Statistics.Classification (ClassificationScore(..))
import Statistics.Classification.ROC (ROC(..), roc)
import System.Directory (doesFileExist)
import System.FilePath ((</>), (<.>))
import qualified Data.ByteString.Lazy as LB
import qualified Data.Vector as V

checkAUC :: V.Vector (ClassificationScore Bool) -> Bool
checkAUC v = 0 <= auc && auc <= 1
  where auc = rocAUC $ roc v

newtype CScore = CScore { unCS :: ClassificationScore Bool }

instance FromNamedRecord CScore where
  parseNamedRecord o = fmap CScore $ ClassificationScore <$> fmap (== ("up" :: String)) (o .: "label") <*> o .: "prob_up"

testExamples :: TestTree
testExamples = testGroup "examples" $ map (\p -> testCase p (ex p)) validp
  where
    validp = filter (/= "ri_czc-2") $ pth <$> inst <*> [1..5 :: Int]
    inst = ["b_dce", "j_dce", "ri_czc"]
    pth prod n = prod ++ "-" ++ show n
    ex b =
      let p = "data" </> b <.> "csv" <.> "xz"
      in doesFileExist p >>= \case
      True -> decodeByName . decompress <$> LB.readFile p >>= \case
        Left err -> assertFailure err
        Right (_, cs) -> assertBool "auc in range" $ checkAUC (V.map unCS cs)
      False -> assertBool "no file" True

propAucRange :: [ClassificationScore Bool] -> Property
propAucRange l = length l >= 1 ==> checkAUC (V.fromList l)

testRoc :: TestTree
testRoc = testGroup "roc" [
    testProperty "auc range" propAucRange
  , testExamples
  ]

max_stars_repo_name: tsbattman/rochs

max_issues_repo_name: tsbattman/rochs

max_forks_repo_name: tsbattman/rochs


text:
{-# LANGUAGE OverloadedStrings #-}
module Numeric.Clustering.Perplexity where

import qualified Data.Massiv.Array as MA
import qualified Numeric.NLOPT as NLOPT
import qualified Numeric.LinearAlgebra as LA
import qualified Control.Monad.Catch as X
import qualified Data.Text as Text
import Numeric (log)

pGivenBeta :: MA.Source r MA.Ix1 Double => MA.Ix1 -> MA.Vector r Double -> Double -> MA.Vector MA.D Double
pGivenBeta rowIndex distances beta =
  let scale ci d = if (ci == rowIndex) then 0 else exp (-d * beta)
      raw = MA.imap scale distances
      sumRaw = MA.sum raw
  in raw MA..* (1 / sumRaw)

safe_xLogx :: Double -> Double
safe_xLogx x = if x > 1e-7 then x * log x else 0

entropy :: MA.Source r MA.Ix1 Double => MA.Vector r Double -> MA.Ix1 -> Double -> Double
entropy distances rowIndex beta = MA.sum $ MA.map (\x -> negate $ safe_xLogx x) $ pGivenBeta rowIndex distances beta

findBeta :: MA.Source r MA.Ix1 Double => Double -> MA.Ix1 -> MA.Vector r Double -> Either Text.Text Double
findBeta perplexity rowIndex distances =
  let targetEntropy = log perplexity
      obj x = abs (targetEntropy - entropy distances rowIndex (x `LA.atIndex` 0))
      algo = NLOPT.NELDERMEAD obj [] Nothing
      stopWhenAny = NLOPT.ObjectiveAbsoluteTolerance 1e-4 NLOPT.:|
                    [NLOPT.MaximumEvaluations 50
                    ]
      problem = NLOPT.LocalProblem 1 stopWhenAny algo
      optimizedE = NLOPT.minimizeLocal problem (LA.fromList [1])
      textError sol = (Text.pack $ show $ NLOPT.solutionResult sol)
                      <> ": x = "
                      <> (Text.pack $ show $ NLOPT.solutionParams sol)
                      <> "; obj(x) = "
                      <> (Text.pack $ show $ NLOPT.solutionCost sol)
  in case optimizedE of
    Left r -> Left $ Text.pack $ show r
    Right s -> case NLOPT.solutionResult s of
      NLOPT.FTOL_REACHED -> Right $ (`LA.atIndex` 0) $ NLOPT.solutionParams s
      _ -> Left $ textError s

data SolveException = SolveException Text.Text deriving (Show)
instance X.Exception SolveException

probabilities :: MA.MonadThrow m => Double -> MA.Matrix MA.U Double -> m (MA.Matrix MA.U Double)
probabilities perplexity distances = do
  let distanceRows = MA.computeAs MA.B $ MA.outerSlices distances
      betaSolsE :: Either Text.Text (MA.Vector MA.U Double)
      betaSolsE = MA.itraverseA (\ix v -> findBeta perplexity ix v) distanceRows
  betas <- case betaSolsE of
    Left err -> X.throwM (SolveException err)
    Right x -> return x
  fmap MA.compute $ MA.stackOuterSlicesM $ MA.izipWith pGivenBeta distanceRows betas

max_stars_repo_name: adamConnerSax/data-science-utils

max_issues_repo_name: adamConnerSax/data-science-utils

max_forks_repo_name: adamConnerSax/data-science-utils