text:
{-# LANGUAGE DataKinds, RankNTypes, TypeFamilies #-}
module TestCommon ((+\), (*\), (**\),
                   listsToParameters,
                   cmpTwo, cmpTwoSimple,
                   qcPropDom, quickCheckTest0, fquad, quad
                  ) where

import Prelude

import qualified Data.Array.DynamicS as OT
import qualified Data.Array.ShapedS as OS
import qualified Data.Vector.Generic as V
import           Numeric.LinearAlgebra (Vector)
import qualified Numeric.LinearAlgebra as HM
import           Test.Tasty
import           Test.Tasty.QuickCheck

import HordeAd hiding (sumElementsVectorOfDual)
import HordeAd.Core.DualClass (Dual)

(+\) :: DualMonad d r m
     => DualNumber d r -> DualNumber d r -> m (DualNumber d r)
(+\) u v = returnLet $ u + v

(*\) :: DualMonad d r m
     => DualNumber d r -> DualNumber d r -> m (DualNumber d r)
(*\) u v = returnLet $ u * v

(**\) :: DualMonad d r m
      => DualNumber d r -> DualNumber d r -> m (DualNumber d r)
(**\) u v = returnLet $ u ** v

-- Checks if 2 numbers are close enough.
close1 :: forall r. (Ord r, Fractional r)
          => r -> r -> Bool
close1 a b = abs (a - b) <= 1e-4

-- Checks if 2 number pairs are close enough.
close2 :: forall r. (Ord r, Fractional r)
          => (r,r) -> (r,r) -> Property
close2 (a1, b1) (a2, b2) = close1 a1 a2 .&&. close1 b1 b2

quad :: DualMonad d r m
     => DualNumber d r -> DualNumber d r -> m (DualNumber d r)
quad x y = do
  x2 <- returnLet $ square x
  y2 <- y *\ y
  tmp <- x2 +\ y2
  tmp +\ 5

fquad :: forall r d m. DualMonad d r m
      => DualNumberVariables d r -> m (DualNumber d r)
fquad variables = do
  let x = var0 variables 0
      y = var0 variables 1
  quad x y

listsToParameters :: forall r. (OT.Storable r)
                  => ([r], [r]) -> Domains r
listsToParameters (a0, a1) =
  (V.fromList a0, V.singleton $ V.fromList a1, V.empty, V.empty)

listsToParameters4 :: ([Double], [Double], [Double], [Double]) -> Domains Double
listsToParameters4 (a0, a1, a2, aX) =
  ( V.fromList a0
  , V.singleton $ V.fromList a1
  , if null a2 then V.empty else V.singleton $ HM.matrix 1 a2
  , if null aX then V.empty else V.singleton $ OT.fromList [length aX] aX )

quickCheckTest0 :: TestName
       -> (forall d r m. ( DualMonad d r m
                         , Floating (Out (DualNumber d (Vector r)))
                         , Floating (Out (DualNumber d (OS.Array '[2] r))) )
           => DualNumberVariables d r -> m (DualNumber d r))
       -> ((Double, Double, Double) -> ([Double], [Double], [Double], [Double]))
       -> TestTree
quickCheckTest0 txt f fArg =
  qcTestRanges txt f (listsToParameters4 . fArg) ((-2, -2, -2), (2, 2, 2)) ((-1e-7, -1e-7, -1e-7), (1e-7, 1e-7, 1e-7)) (-10, 10)

-- A quick check to compare the derivatives and values of 2 given functions.
cmpTwo
  :: ( m ~ DualMonadForward r, d ~ 'DModeDerivative, Dual d r ~ r
     , DualMonad d r m )
  => (DualNumberVariables d r -> m (DualNumber d r))
  -> (DualNumberVariables d r -> m (DualNumber d r))
  -> Domains r
  -> Domains r
  -> Domains r
  -> Domains r
  -> Property
cmpTwo f1 f2 params1 params2 ds1 ds2 =
  close2 (dFastForward f1 params1 ds1) (dFastForward f2 params2 ds2)

-- A quick check to compare the derivatives and values of 2 given functions.
cmpTwoSimple
  :: ( m ~ DualMonadForward r, d ~ 'DModeDerivative, Dual d r ~ r
     , DualMonad d r m )
  => (DualNumberVariables d r -> m (DualNumber d r))
  -> (DualNumberVariables d r -> m (DualNumber d r))
  -> Domains r
  -> Domains r
  -> Property
cmpTwoSimple f1 f2 parameters ds =
  cmpTwo f1 f2 parameters parameters ds ds

-- A quick consistency check of all the kinds of derivatives and gradients
-- and all kinds of computing the value of the objective function.
qcPropDom :: (forall d r m. ( DualMonad d r m
                            , r ~ Double
                            , Floating (Out (DualNumber d (Vector r)))
                            , Floating (Out (DualNumber d (OS.Array '[2] r))) )
              => DualNumberVariables d r -> m (DualNumber d r))
       -> Domains Double
       -> Domains Double
       -> Domains Double
       -> Double
       -> Property
qcPropDom f args ds perturbation dt =
      let ff@(derivative, ffValue) = dFastForward f args ds
          (derivativeAtPerturbation, valueAtPerturbation) = dFastForward f args perturbation
          (gradient, revValue) = dReverse dt f args
      in -- Two forward derivative implementations agree fully:
         dForward f args ds === ff
         -- Objective function value from gradients is the same.
         .&&. ffValue == revValue
         -- Gradients and derivatives agree.
         .&&. close1 (dt * derivative)
                     (dotParameters gradient ds)
         -- Objective function value is unaffected by perturbation.
         .&&. ffValue == valueAtPerturbation
         -- Derivative approximates the perturbation of value.
         .&&. close1 (primalValue
                                  f (addParameters
                                                   args perturbation))
                     (ffValue + derivativeAtPerturbation)

-- A quick consistency check of all the kinds of derivatives and gradients
-- and all kinds of computing the value of the objective function.
qcPropFArg :: (forall d r m. ( DualMonad d r m
                         , Floating (Out (DualNumber d (Vector r)))
                         , Floating (Out (DualNumber d (OS.Array '[2] r))) )
               => DualNumberVariables d r -> m (DualNumber d r))
       -> ((Double, Double, Double) -> Domains Double)
       -> (Double, Double, Double)
       -> (Double, Double, Double)
       -> (Double, Double, Double)
       -> Double
       -> Property
qcPropFArg f fArgDom xyz dsRaw perturbationRaw dt =
      let args = fArgDom xyz
          ds = fArgDom dsRaw
          perturbation = fArgDom perturbationRaw
      in qcPropDom f args ds perturbation dt

-- A quick consistency check of all the kinds of derivatives and gradients
-- and all kinds of computing the value of the objective function.
qcTestRanges :: TestName
       -> (forall d r m. ( DualMonad d r m
                         , Floating (Out (DualNumber d (Vector r)))
                         , Floating (Out (DualNumber d (OS.Array '[2] r))) )
           => DualNumberVariables d r -> m (DualNumber d r))
       -> ((Double, Double, Double) -> Domains Double)
       -> ((Double, Double, Double), (Double, Double, Double))
       -> ((Double, Double, Double), (Double, Double, Double))
       -> (Double, Double)
       -> TestTree
qcTestRanges txt f fArgDom dsRange perturbationRange dtRange =
  testProperty txt $
  forAll (choose dsRange) $ \xyz dsRaw ->
  forAll (choose perturbationRange) $ \perturbationRaw ->
  forAll (choose dtRange) $ \dt ->
  qcPropFArg f fArgDom xyz dsRaw perturbationRaw dt

max_stars_repo_name: Mikolaj/horde-ad

max_issues_repo_name: Mikolaj/horde-ad

max_forks_repo_name: Mikolaj/horde-ad


text:
--
-- Spec file for luSolve.hs
--
-- Tests the exposed LU decomposition and linear system solver.
--
--
module Numeric.LinearAlgebra.LUSolveSpec where

import           Test.Hspec

import           Data.Matrix.Generic           as M hiding (zipWith)
import qualified Data.Vector.Unboxed           as V
import           Numeric.LinearAlgebra.LUSolve
import           System.Random

-- A useful, but strangely non-standard function: given an integer n
-- and a list, break the input list up into lists of length n and
-- return the list of those lists.
--
-- The name reflects that it "bundles" a list into lists of length n.
--
-- If the length of the input list is not a multiple of n, the final
-- sublist will be shorter than n elements.
--
-- The definition here is works for infinite lists, which is perhaps
-- the most common use-case for this function.
--
bundle :: Int -> [ a ] -> [[ a ]]
bundle _ [] = []
bundle n xs = take n xs : bundle n (drop n xs)

mVals :: [ Double ]
mVals = randoms (mkStdGen 1)

-- Make an infinite list of random n * n matrices.
--
randomSquareMatrices :: Int -> [ M.Matrix V.Vector Double ]
randomSquareMatrices n = Prelude.map (\vs -> M.fromLists (bundle n vs)) (bundle (n * n) mVals)

vVals :: [ Double ]
vVals = randoms (mkStdGen 2)

-- Make an infinite list of random n-vectors.
--
randomColumnVectors :: Int -> [ M.Matrix V.Vector Double ]
randomColumnVectors n = Prelude.map (\vs -> M.fromLists (bundle 1 vs )) (bundle n vVals)

-- Solve the linear system Ax = b for x, then compute the maximum
-- absolute value of the vector of differences Ax - b.  The later
-- is the value of the function.
--
checkLinearSystemSolution :: M.Matrix V.Vector Double  -- matrix A
                          -> M.Matrix V.Vector Double  -- right hand vector b
                          -> Double                    -- max (abs (Ax - b))
checkLinearSystemSolution a b = let
    x  = luSolve (luFactor a) b
    b' = matrixMultiply a x
    in
      maxAbsDiff b b'

-- Multiply two matrices
--
matrixMultiply :: M.Matrix V.Vector Double
               -> M.Matrix V.Vector Double
               -> M.Matrix V.Vector Double
matrixMultiply a b = let
    (ra, ca) = dim a
    (rb, cb) = dim b
    in
      if ca /= rb
      then error "incompatible dimensions in matrixMultiply"
      else fromLists [[ sum [ a ! (i, k) * b ! (k, j) | k <- [0 .. (ca - 1)]]
                                                      | j <- [0 .. (cb - 1)]]
                                                      | i <- [0 .. (ra - 1)]]

-- Return the largest difference (by absolute value) between
-- two matrices.
--
maxAbsDiff :: M.Matrix V.Vector Double
           -> M.Matrix V.Vector Double
           -> Double
maxAbsDiff (M.Matrix m n _ _ v) (M.Matrix m' n' _ _ v') =
    if m == m' && n == n'
    then V.maximum $ V.zipWith (\y y' -> abs (y - y')) v v'
    else error "unequal matrix dimensions in maxAbsDiff"

-- Return the largest difference (by relative value) between
-- two matrices.
--
maxRelDiff :: M.Matrix V.Vector Double
           -> M.Matrix V.Vector Double
           -> Double
maxRelDiff (Matrix m n _ _ v) (Matrix m' n' _ _ v') =
    if m == m' && n == n'
    then V.maximum $ V.zipWith (\y y' -> abs (2 * (y - y') / (abs y + abs y'))) v v'
    else error "unequal matrix dimensions in maxRelDiff"

-- Threshold for accepting max (abs (Ax - b)) as correct.
--
matrixSpecAbsDiff :: Double
matrixSpecAbsDiff = 1.0e-9

-- Generate n linear systems of size 'size' and solve them by luSolve.
-- Return a list of the maximum absolute value of the error (max (abs (Ax - b)))
-- for each case.
--
runCheck :: Int         -- number of cases to run
         -> Int         -- number of equations in each linear system
         -> [ Double ]  -- list of maximum absolute values of the errors
runCheck n size = take n $ zipWith checkLinearSystemSolution (randomSquareMatrices size)
                                                             (randomColumnVectors  size)

-- The tests themselves
--
spec :: Spec
spec = do
    -- describe "Numeric.LinearAlgebra.LUSolve.luSolve" $ do
    --     it "Solve 1000 random 10 * 10 systems" $ do
    --         (maximum $ runCheck 1000 10)
    --             `shouldSatisfy` (< matrixSpecAbsDiff)

    -- describe "Numeric.LinearAlgebra.LUSolve.luSolve" $ do
    --     it "Solve 100 random 50 * 50 systems" $ do
    --         (maximum $ runCheck 100 50)
    --             `shouldSatisfy` (< matrixSpecAbsDiff)

    -- describe "Numeric.LinearAlgebra.LUSolve.luSolve" $ do
    --     it "Solve 100 random 100 * 100 systems" $ do
    --         (maximum $ runCheck 100 100)
    --             `shouldSatisfy` (< matrixSpecAbsDiff)

    -- describe "Numeric.LinearAlgebra.LUSolve.luSolve" $ do
    --     it "Solve 100 random 500 * 500 systems" $ do
    --         (maximum $ runCheck 100 500)
    --             `shouldSatisfy` (< matrixSpecAbsDiff)

    describe "Numeric.LinearAlgebra.LUSolve.luSolve" $ do
        it "Solve 100 random 1000 * 1000 systems" $ do
            (maximum $ runCheck 100 1000)
                `shouldSatisfy` (< matrixSpecAbsDiff)

main :: IO ()
main = hspec spec

ext: hs


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE Strict #-}
{-# LANGUAGE FlexibleContexts #-}

module CPVO.IO.Fortran (
  loadFile
  ) where

--import Numeric.LinearAlgebra
--import Numeric.LinearAlgebra.Data hiding (find)
--import Data.List (findIndex,groupBy)
--import Data.Maybe (fromJust)
--import qualified Language.C.Inline as C

loadFile :: IO ()
loadFile = do
  putStrLn "=======done : loadFile@CPVO.IO.Fortran"

max_stars_repo_name: hasanalrasyid/cpvoh

max_issues_repo_name: hasanalrasyid/cpvoh

max_forks_repo_name: hasanalrasyid/cpvoh


text:
{-# LANGUAGE DeriveAnyClass #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ImportQualifiedPost #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE ViewPatterns #-}
{-# OPTIONS_GHC -Wwarn #-}

module Data.CDF
  ( Centile(..)
  , renderCentile
  , briefCentiles
  , stdCentiles
  , nEquicentiles
  , CDFError (..)
  , CDF(..)
  , CDFIx (..)
  , KnownCDF (..)
  , liftCDFVal
  , unliftCDFVal
  , centilesCDF
  , filterCDF
  , zeroCDF
  , projectCDF
  , projectCDF'
  , indexCDF
  , DirectCDF
  , cdf
  , mapToCDF
  , Divisible (..)
  , Combine (..)
  , stdCombine1
  , stdCombine2
  , CDF2
  , collapseCDFs
  , cdf2OfCDFs
  --
  , module Data.SOP.Strict
  ) where

import Prelude (String, (!!), error, head, show)
import Cardano.Prelude hiding (head, show)

import Data.Aeson (FromJSON(..), ToJSON(..))
import Data.SOP.Strict
import Data.Time.Clock (NominalDiffTime)
import Data.Vector qualified as Vec
import Statistics.Sample qualified as Stat
import Text.Printf (printf)

import Ouroboros.Consensus.Util.Time (secondsToNominalDiffTime)

-- | Centile specifier: a fractional in range of [0; 1].
newtype Centile =
  Centile { unCentile :: Double }
  deriving (Eq, Generic, FromJSON, ToJSON, Show)
  deriving anyclass NFData

renderCentile :: Int -> Centile -> String
renderCentile width = \case
  Centile x    -> printf ("%0."<>show (width-2)<>"f") x

briefCentiles :: [Centile]
briefCentiles =
  [ Centile 0.5, Centile 0.9, Centile 1.0 ]

stdCentiles :: [Centile]
stdCentiles =
  [ Centile 0.01, Centile 0.05
  , Centile 0.1, Centile 0.2, Centile 0.3, Centile 0.4
  , Centile 0.5, Centile 0.6
  , Centile 0.7, Centile 0.75
  , Centile 0.8, Centile 0.85, Centile 0.875
  , Centile 0.9, Centile 0.925, Centile 0.95, Centile 0.97, Centile 0.98, Centile 0.99
  , Centile 0.995, Centile 0.997, Centile 0.998, Centile 0.999
  , Centile 0.9995, Centile 0.9997, Centile 0.9998, Centile 0.9999
  ]

-- | Given a N-large population, produce centiles for each element, except for min and max.
--   We don't need min and max, because CDF already has range.
nEquicentiles :: Int -> [Centile]
nEquicentiles n =
  if reindices == indices
  then Centile <$> centiles
  else error $ printf "centilesForN:  reindices for %d: %s, indices: %s" n (show reindices) (show indices)
 where
   reindices = centiles <&> runCentile n
   centiles  = [ step * (fromIntegral i + 0.5) | i <- indices ]
   indices   = if n > 2
               then [1 .. n - 2] -- ignore first and last indices, standing for min and max.
               else []
   step :: Double
   step = 1.0 / fromIntegral n

-- | Given a centile of N-large population, produce index of the population element referred by centile.
{-# INLINE runCentile #-}
runCentile :: Int -> Double -> Int
runCentile n centile = floor (fromIntegral n * centile)
                       & min (n - 1)

{-# INLINE vecCentile #-}
vecCentile :: Vec.Vector a -> Int -> Centile -> a
vecCentile vec n (Centile c) = vec Vec.! runCentile n c

--
-- * Parametric CDF (cumulative distribution function)
--
data CDF p a =
  CDF
  { cdfSize      :: Int
  , cdfAverage   :: Double
  , cdfStddev    :: Double
  , cdfRange     :: (a, a)
  , cdfSamples  :: [(Centile, p a)]
  }
  deriving (Eq, Functor, Generic, Show)
  deriving anyclass NFData

instance (FromJSON (p a), FromJSON a) => FromJSON (CDF p a)
instance (  ToJSON (p a),   ToJSON a) => ToJSON   (CDF p a)

-- * Singletons
--
data CDFIx p where
  CDFI :: CDFIx I
  CDF2 :: CDFIx (CDF I)

class KnownCDF a where
  cdfIx :: CDFIx a

instance KnownCDF      I  where cdfIx = CDFI
instance KnownCDF (CDF I) where cdfIx = CDF2

centilesCDF :: CDF p a -> [Centile]
centilesCDF = fmap fst . cdfSamples

zeroCDF :: (Num a) => CDF p a
zeroCDF =
  CDF
  { cdfSize    = 0
  , cdfAverage = 0
  , cdfStddev  = 0
  , cdfRange   = (0, 0)
  , cdfSamples = mempty
  }

filterCDF :: ((Centile, p a) -> Bool) -> CDF p a -> CDF p a
filterCDF f d =
  d { cdfSamples = cdfSamples d & filter f }

indexCDF :: Int -> CDF p a -> p a
indexCDF i d = snd $ cdfSamples d !! i

projectCDF :: Centile -> CDF p a -> Maybe (p a)
projectCDF p = fmap snd . find ((== p) . fst) . cdfSamples

projectCDF' :: String -> Centile -> CDF p a -> p a
projectCDF' desc p =
  maybe (error er) snd . find ((== p) . fst) . cdfSamples
 where
   er = printf "Missing centile %f in %s" (show $ unCentile p) desc

--
-- * Trivial instantiation: samples are value points
--
type DirectCDF a = CDF I a

liftCDFVal :: Real a => a -> CDFIx p -> p a
liftCDFVal x = \case
  CDFI -> I x
  CDF2 -> CDF { cdfSize    = 1
              , cdfAverage = fromRational (toRational x)
              , cdfStddev  = 0
              , cdfRange   = (x, x)
              , cdfSamples = []
              , .. }

unliftCDFVal :: Divisible a => CDFIx p -> p a -> a
unliftCDFVal CDFI (I x) =x
unliftCDFVal CDF2 CDF{cdfRange = (mi, ma)} = (mi + ma) `divide` 2

cdf :: (Real a, KnownCDF p) => [Centile] -> [a] -> CDF p a
cdf centiles (sort -> sorted) =
  CDF
  { cdfSize        = size
  , cdfAverage     = Stat.mean   doubleVec
  , cdfStddev      = Stat.stdDev doubleVec
  , cdfRange       = (mini, maxi)
  , cdfSamples =
    (    (Centile   0, liftCDFVal mini ix) :) .
    (<> [(Centile 1.0, liftCDFVal maxi ix) ]) $
    centiles <&>
      \spec ->
        let sample = if size == 0 then 0
                     else vecCentile vec size spec
        in (,) spec (liftCDFVal sample ix)
  }
 where ix          = cdfIx
       vec         = Vec.fromList sorted
       size        = length vec
       doubleVec   = fromRational . toRational <$> vec
       (,) mini maxi =
         if size == 0
         then (0,           0)
         else (vec Vec.! 0, Vec.last vec)

mapToCDF :: Real a => (b -> a) -> [Centile] -> [b] -> DirectCDF a
mapToCDF f pspecs xs = cdf pspecs (f <$> xs)

type CDF2 a = CDF (CDF I) a

data CDFError
  = CDFIncoherentSamplingLengths  [Int]
  | CDFIncoherentSamplingCentiles [[Centile]]
  | CDFEmptyDataset
  deriving Show

-- | Avoiding `Fractional`
class Real a => Divisible a where
  divide :: a -> Double -> a

instance Divisible Double where
  divide = (/)

instance Divisible Int where
  divide x by = round $ fromIntegral x / by

instance Divisible Integer where
  divide x by = round $ fromIntegral x / by

instance Divisible Word32 where
  divide x by = round $ fromIntegral x / by

instance Divisible Word64 where
  divide x by = round $ fromIntegral x / by

instance Divisible NominalDiffTime where
  divide x by = x / secondsToNominalDiffTime by

-- * Combining population stats
data Combine p a
  = Combine
    { cWeightedAverages :: !([(Int, Double)] -> Double)
    , cStddevs          :: !([Double] -> Double)
    , cRanges           :: !([(a, a)] -> (a, a))
    , cWeightedSamples  :: !([(Int, a)] -> a)
    , cCDF              :: !([p a] -> Either CDFError (CDF I a))
    }

stdCombine1 :: forall a. (Divisible a) => [Centile] -> Combine I a
stdCombine1 cs =
  Combine
  { cWeightedAverages = weightedAverage
  , cRanges           = outerRange
  , cStddevs          = maximum          -- it's an approximation
  , cWeightedSamples  = weightedAverage
  , cCDF              = Right . cdf cs . fmap unI
  }
  where
    weightedAverage :: forall b. (Divisible b) => [(Int, b)] -> b
    weightedAverage xs = (`divide` (fromIntegral . sum $ fst <$> xs)) . sum $
                          xs <&> \(size, avg) -> fromIntegral size * avg
    outerRange      xs = (,) (minimum $ fst <$> xs)
                             (maximum $ snd <$> xs)

stdCombine2 :: Divisible a => [Centile] -> Combine (CDF I) a
stdCombine2 cs =
  let c@Combine{..} = stdCombine1 cs in
  Combine
  { cCDF = collapseCDFs c
  , ..
  }

-- | Collapse:  Given a ([Value] -> CDF I) function, and a list of (CDF I), produce a (CDF I).
--
collapseCDFs :: forall a. Combine I a -> [CDF I a] -> Either CDFError (CDF I a)
collapseCDFs _ [] = Left CDFEmptyDataset
collapseCDFs Combine{..} xs = do
  unless (all (head lengths ==) lengths) $
    Left $ CDFIncoherentSamplingLengths lengths
  unless (null incoherent) $
    Left $ CDFIncoherentSamplingCentiles (fmap fst <$> incoherent)
  pure CDF
    { cdfSize    = sum sizes
    , cdfAverage = xs <&> cdfAverage & cWeightedAverages . zip sizes
    , cdfRange   = xs <&> cdfRange   & cRanges
    , cdfStddev  = xs <&> cdfStddev  & cStddevs
    , cdfSamples = coherent <&>
                   fmap (I . cWeightedSamples . zip sizes . fmap unI)
    }
 where
   sizes   = xs <&> cdfSize
   samples = xs <&> cdfSamples
   lengths = length <$> samples

   centileOrdered :: [[(Centile, I a)]] -- Each sublist must (checked) have the same Centile.
   centileOrdered = transpose samples

   coherent :: [(Centile, [I a])]
   (incoherent, coherent) = partitionEithers $ centileOrdered <&>
     \case
       [] -> error "cdfOfCDFs:  empty list of centiles, hands down."
       xxs@((c, _):(fmap fst -> cs)) -> if any (/= c) cs
                                        then Left xxs
                                        else Right (c, snd <$> xxs)

-- | Polymorphic, but practically speaking, intended for either:
--    1. given a ([I]     -> CDF I) function, and a list of (CDF I),       produce a CDF (CDF I), or
--    2. given a ([CDF I] -> CDF I) function, and a list of (CDF (CDF I)), produce a CDF (CDF I)
cdf2OfCDFs :: forall a p. Combine p a -> [CDF p a] -> Either CDFError (CDF (CDF I) a)
cdf2OfCDFs _ [] = Left CDFEmptyDataset
cdf2OfCDFs Combine{..} xs = do
  unless (all (head lengths ==) lengths) $
    Left $ CDFIncoherentSamplingLengths lengths
  unless (null incoherent) $
    Left $ CDFIncoherentSamplingCentiles (fmap fst <$> incoherent)

  cdfSamples <- mapM sequence -- ..to  Either CDFError [(Centile, CDF I a)]
                  (coherent <&> fmap cCDF :: [(Centile, Either CDFError (CDF I a))])
  pure CDF
    { cdfSize    = sum sizes
    , cdfAverage = xs <&> cdfAverage & cWeightedAverages . zip sizes
    , cdfRange   = xs <&> cdfRange   & cRanges
    , cdfStddev  = xs <&> cdfStddev  & cStddevs
    , cdfSamples = cdfSamples
    }
 where
   sizes   = xs <&> cdfSize
   samples = xs <&> cdfSamples
   lengths = length <$> samples

   centileOrdered :: [[(Centile, p a)]]
   centileOrdered = transpose samples

   (incoherent, coherent) = partitionEithers $ centileOrdered <&>
     \case
       [] -> error "cdfOfCDFs:  empty list of centiles, hands down."
       xxs@((c, _):(fmap fst -> cs)) -> if any (/= c) cs
                                        then Left xxs
                                        else Right (c, snd <$> xxs)


text:
-- |
-- Module      : Southpaw.Interactive.Lenses
-- Description :
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created October 15 2015

-- TODO | -
--        -

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC Pragmas
--------------------------------------------------------------------------------------------------------------------------------------------

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module Southpaw.Interactive.Lenses where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Data.IORef
import Data.Complex
import Data.Functor ((<$>))
import Control.Lens

import Graphics.UI.Gtk
import qualified Graphics.Rendering.Cairo as Cairo

import qualified Southpaw.Interactive.Types as App

--------------------------------------------------------------------------------------------------------------------------------------------
-- Lenses
--------------------------------------------------------------------------------------------------------------------------------------------

-- App -------------------------------------------------------------------------------------------------------------------------------------

window :: Lens (App.App layout state) (App.App layout state) Window Window
window f s = (\new -> s { App._window=new }) <$> f (App._window s)

layout :: Lens (App.App layout state) (App.App layout state) layout layout
layout f s = (\new -> s { App._layout=new }) <$> f (App._layout s)

-- TODO: Make this a getter (since setting this property won't affect the windw size) (?)
size :: Lens (App.App layout state) (App.App layout state) (Complex Int) (Complex Int)
size f s = (\new -> s { App._appsize=new }) <$> f (App._appsize s)

state :: Lens (App.App layout state) (App.App layout state) (IORef state) (IORef state)
state f s = (\new -> s { App._state=new }) <$> f (App._state s)

-- CanvasLayout ----------------------------------------------------------------------------------------------------------------------------

canvas :: Lens (App.FullCanvasLayout) (App.FullCanvasLayout) (DrawingArea) (DrawingArea)
canvas f s = (\new -> s { App._canvas=new }) <$> f (App._canvas s)

-- Settings --------------------------------------------------------------------------------------------------------------------------------

-- Generic ---------------------------------------------------------------------------------------------------------------------------------

max_stars_repo_name: SwiftsNamesake/Southpaw

max_issues_repo_name: SwiftsNamesake/Southpaw

max_forks_repo_name: SwiftsNamesake/Southpaw


text:
module Main where

import Options.Applicative

import Pipes
import qualified Pipes.Prelude as P

import System.Random.MWC

import Statistics.Quantile.Bench.Gen (sampleGamma)

data Params = Params Int Double Double
  deriving (Eq, Show)

params = Params <$> argument auto (metavar "COUNT")
                <*> argument auto (metavar "ALPHA")
                <*> argument auto (metavar "BETA")

args = info (helper <*> params)
            (   fullDesc
             <> progDesc "Write COUNT real numbers drawn from a gamma distribution."
            )

main :: IO ()
main = do
    ps <- execParser args
    gen <- createSystemRandom
    runEffect $ gamma gen ps >-> P.map show >-> P.stdoutLn
  where
    gamma gen (Params n alpha beta) = P.replicateM n $ sampleGamma gen alpha beta

max_stars_repo_name: fractalcat/slides

max_issues_repo_name: fractalcat/slides

max_forks_repo_name: fractalcat/slides


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE DefaultSignatures #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE FunctionalDependencies #-}
#ifdef TRUSTWORTHY
{-# LANGUAGE Trustworthy #-}
#endif

#include "lens-common.h"

-----------------------------------------------------------------------------
-- |
-- Module      :  Control.Lens.Each
-- Copyright   :  (C) 2012-16 Edward Kmett
-- License     :  BSD-style (see the file LICENSE)
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  experimental
-- Portability :  non-portable
--
-----------------------------------------------------------------------------
module Control.Lens.Each
  (
  -- * Each
    Each(..)
  ) where

import Prelude ()

import Control.Lens.Traversal
import Control.Lens.Internal.ByteString
import Control.Lens.Internal.Prelude
import Data.Array.Unboxed as Unboxed
import Data.Array.IArray as IArray
import Data.ByteString as StrictB
import Data.ByteString.Lazy as LazyB
import Data.Complex
import Data.HashMap.Lazy as HashMap
import Data.IntMap as IntMap
import Data.Map as Map
import Data.Sequence as Seq
import Data.Text.Lens (text)
import Data.Text as StrictT
import Data.Text.Lazy as LazyT
import Data.Tree as Tree
import Data.Vector.Generic.Lens (vectorTraverse)
import qualified Data.Vector as Vector
import qualified Data.Vector.Primitive as Prim
import Data.Vector.Primitive (Prim)
import qualified Data.Vector.Storable as Storable
import Data.Vector.Storable (Storable)
import qualified Data.Vector.Unboxed as Unboxed
import Data.Vector.Unboxed (Unbox)
import Data.Word

-- $setup
-- >>> :set -XNoOverloadedStrings
-- >>> import Control.Lens
-- >>> import Data.Text.Strict.Lens as Text
-- >>> import Data.Char as Char

-- | Extract 'each' element of a (potentially monomorphic) container.
--
-- Notably, when applied to a tuple, this generalizes 'Control.Lens.Traversal.both' to arbitrary homogeneous tuples.
--
-- >>> (1,2,3) & each *~ 10
-- (10,20,30)
--
-- It can also be used on monomorphic containers like 'StrictT.Text' or 'StrictB.ByteString'.
--
-- >>> over each Char.toUpper ("hello"^.Text.packed)
-- "HELLO"
--
-- >>> ("hello","world") & each.each %~ Char.toUpper
-- ("HELLO","WORLD")
class Each s t a b | s -> a, t -> b, s b -> t, t a -> s where
  each :: Traversal s t a b
  default each :: (Traversable g, s ~ g a, t ~ g b) => Traversal s t a b
  each = traverse
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' (a,a) (b,b) a b@
instance (a~a', b~b') => Each (a,a') (b,b') a b where
  each f ~(a,b) = (,) <$> f a <*> f b
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' (a,a,a) (b,b,b) a b@
instance (a~a2, a~a3, b~b2, b~b3) => Each (a,a2,a3) (b,b2,b3) a b where
  each f ~(a,b,c) = (,,) <$> f a <*> f b <*> f c
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' (a,a,a,a) (b,b,b,b) a b@
instance (a~a2, a~a3, a~a4, b~b2, b~b3, b~b4) => Each (a,a2,a3,a4) (b,b2,b3,b4) a b where
  each f ~(a,b,c,d) = (,,,) <$> f a <*> f b <*> f c <*> f d
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' (a,a,a,a,a) (b,b,b,b,b) a b@
instance (a~a2, a~a3, a~a4, a~a5, b~b2, b~b3, b~b4, b~b5) => Each (a,a2,a3,a4,a5) (b,b2,b3,b4,b5) a b where
  each f ~(a,b,c,d,e) = (,,,,) <$> f a <*> f b <*> f c <*> f d <*> f e
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' (a,a,a,a,a,a) (b,b,b,b,b,b) a b@
instance (a~a2, a~a3, a~a4, a~a5, a~a6, b~b2, b~b3, b~b4, b~b5, b~b6) => Each (a,a2,a3,a4,a5,a6) (b,b2,b3,b4,b5,b6) a b where
  each f ~(a,b,c,d,e,g) = (,,,,,) <$> f a <*> f b <*> f c <*> f d <*> f e <*> f g
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' (a,a,a,a,a,a,a) (b,b,b,b,b,b,b) a b@
instance (a~a2, a~a3, a~a4, a~a5, a~a6, a~a7, b~b2, b~b3, b~b4, b~b5, b~b6, b~b7) => Each (a,a2,a3,a4,a5,a6,a7) (b,b2,b3,b4,b5,b6,b7) a b where
  each f ~(a,b,c,d,e,g,h) = (,,,,,,) <$> f a <*> f b <*> f c <*> f d <*> f e <*> f g <*> f h
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' (a,a,a,a,a,a,a,a) (b,b,b,b,b,b,b,b) a b@
instance (a~a2, a~a3, a~a4, a~a5, a~a6, a~a7, a~a8, b~b2, b~b3, b~b4, b~b5, b~b6, b~b7, b~b8) => Each (a,a2,a3,a4,a5,a6,a7,a8) (b,b2,b3,b4,b5,b6,b7,b8) a b where
  each f ~(a,b,c,d,e,g,h,i) = (,,,,,,,) <$> f a <*> f b <*> f c <*> f d <*> f e <*> f g <*> f h <*> f i
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' (a,a,a,a,a,a,a,a,a) (b,b,b,b,b,b,b,b,b) a b@
instance (a~a2, a~a3, a~a4, a~a5, a~a6, a~a7, a~a8, a~a9, b~b2, b~b3, b~b4, b~b5, b~b6, b~b7, b~b8, b~b9) => Each (a,a2,a3,a4,a5,a6,a7,a8,a9) (b,b2,b3,b4,b5,b6,b7,b8,b9) a b where
  each f ~(a,b,c,d,e,g,h,i,j) = (,,,,,,,,) <$> f a <*> f b <*> f c <*> f d <*> f e <*> f g <*> f h <*> f i <*> f j
  {-# INLINE each #-}

-- | @'each' :: ('RealFloat' a, 'RealFloat' b) => 'Traversal' ('Complex' a) ('Complex' b) a b@
instance Each (Complex a) (Complex b) a b where
  each f (a :+ b) = (:+) <$> f a <*> f b
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' ('Map' c a) ('Map' c b) a b@
instance (c ~ d) => Each (Map c a) (Map d b) a b where
  each = traversed
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' ('Map' c a) ('Map' c b) a b@
instance Each (IntMap a) (IntMap b) a b where
  each = traversed
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' ('HashMap' c a) ('HashMap' c b) a b@
instance (c ~ d) => Each (HashMap c a) (HashMap d b) a b where
  each = traversed
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' [a] [b] a b@
instance Each [a] [b] a b where
  each = traversed
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' (NonEmpty a) (NonEmpty b) a b@
instance Each (NonEmpty a) (NonEmpty b) a b

-- | @'each' :: 'Traversal' ('Identity' a) ('Identity' b) a b@
instance Each (Identity a) (Identity b) a b

-- | @'each' :: 'Traversal' ('Maybe' a) ('Maybe' b) a b@
instance Each (Maybe a) (Maybe b) a b

-- | @'each' :: 'Traversal' ('Either' a a) ('Either' b b) a b@
--
-- @since 4.18
instance (a~a', b~b') => Each (Either a a') (Either b b') a b where
  each f (Left a)   = Left <$> f a
  each f (Right a ) = Right <$> f a
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' ('Seq' a) ('Seq' b) a b@
instance Each (Seq a) (Seq b) a b where
  each = traversed
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' ('Tree' a) ('Tree' b) a b@
instance Each (Tree a) (Tree b) a b

-- | @'each' :: 'Traversal' ('Vector.Vector' a) ('Vector.Vector' b) a b@
instance Each (Vector.Vector a) (Vector.Vector b) a b where
  each = vectorTraverse
  {-# INLINE each #-}

-- | @'each' :: ('Prim' a, 'Prim' b) => 'Traversal' ('Prim.Vector' a) ('Prim.Vector' b) a b@
instance (Prim a, Prim b) => Each (Prim.Vector a) (Prim.Vector b) a b where
  each = vectorTraverse
  {-# INLINE each #-}

-- | @'each' :: ('Storable' a, 'Storable' b) => 'Traversal' ('Storable.Vector' a) ('Storable.Vector' b) a b@
instance (Storable a, Storable b) => Each (Storable.Vector a) (Storable.Vector b) a b where
  each = vectorTraverse
  {-# INLINE each #-}

-- | @'each' :: ('Unbox' a, 'Unbox' b) => 'Traversal' ('Unboxed.Vector' a) ('Unboxed.Vector' b) a b@
instance (Unbox a, Unbox b) => Each (Unboxed.Vector a) (Unboxed.Vector b) a b where
  each = vectorTraverse
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' 'StrictT.Text' 'StrictT.Text' 'Char' 'Char'@
instance (a ~ Char, b ~ Char) => Each StrictT.Text StrictT.Text a b where
  each = text
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' 'LazyT.Text' 'LazyT.Text' 'Char' 'Char'@
instance (a ~ Char, b ~ Char) => Each LazyT.Text LazyT.Text a b where
  each = text
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' 'StrictB.ByteString' 'StrictB.ByteString' 'Word8' 'Word8'@
instance (a ~ Word8, b ~ Word8) => Each StrictB.ByteString StrictB.ByteString a b where
  each = traversedStrictTree
  {-# INLINE each #-}

-- | @'each' :: 'Traversal' 'LazyB.ByteString' 'LazyB.ByteString' 'Word8' 'Word8'@
instance (a ~ Word8, b ~ Word8) => Each LazyB.ByteString LazyB.ByteString a b where
  each = traversedLazy
  {-# INLINE each #-}

-- | @'each' :: 'Ix' i => 'Traversal' ('Array' i a) ('Array' i b) a b@
instance (Ix i, i ~ j) => Each (Array i a) (Array j b) a b where
  each f arr = array (bounds arr) <$> traverse (\(i,a) -> (,) i <$> f a) (IArray.assocs arr)
  {-# INLINE each #-}

-- | @'each' :: ('Ix' i, 'IArray' 'UArray' a, 'IArray' 'UArray' b) => 'Traversal' ('Array' i a) ('Array' i b) a b@
instance (Ix i, IArray UArray a, IArray UArray b, i ~ j) => Each (UArray i a) (UArray j b) a b where
  each f arr = array (bounds arr) <$> traverse (\(i,a) -> (,) i <$> f a) (IArray.assocs arr)
  {-# INLINE each #-}

max_stars_repo_name: felixonmars/lens

max_issues_repo_name: fdelacruz-asapp/lens

max_forks_repo_name: fdelacruz-asapp/lens


text:
module Main where
import qualified Data.ByteString.Lazy as BS hiding (putStrLn)
import qualified Data.ByteString.Lazy as BS (putStrLn)
import Codec.Compression.GZip (decompress)
import Numeric.LinearAlgebra hiding (split)
import Prelude hiding (readFile)
import System.Random ( StdGen, mkStdGen )
import System.Random.Shuffle
import Control.Monad
import Data.List.Split
import Data.List

import Control.Monad.State
import Control.Monad (join)
import Data.Bifunctor (bimap)
import Data.ByteString.Conversion.From
import Data.ByteString.Char8 (readInt)
import Data.Int

data Network = Network { layers  :: [Int]
                       , biases  :: [Vector R]
                       , weights :: [Matrix R]
                       , rate    :: Double
                       } deriving (Show)

main :: IO ()
main = do
  let nInputs  = 784
  let nHidden  = 30
  let nOutputs = 10
  let seed = 1
  let layers = [nInputs, nHidden, nOutputs]
  let biases = fmap (\n -> randomVector seed Gaussian n) $ drop 1 layers
  let weights = fmap (\(x, y) -> randomMatrix x y seed + 10)  $ zip (drop 1 layers) layers
  let rate = 1.5
  let net = Network layers biases weights rate

  trainingData <- readImages "images.gz" "labels.gz"
  let vec = snd $ trainingData !! 0
  --print . show . feedForward net $ vec
  print . show . fst $ trainingData !! 1
  --let input = (take 20 trainingData)
  --t <- forM input $ \(digit, input) ->
      --return (digit, train net 20 g input)

------------------- Reading MNIST -------------------------

readImages :: String -> String -> IO ([(Vector R, Vector R)])
readImages iPath lPath = do
  iBytes <- decompress <$> BS.readFile iPath
  lBytes <- decompress <$> BS.readFile lPath
  let is = readImage iBytes <$> [1..50000]
  let ls = readLabel lBytes <$> [1..50000]
  return $ zip ls is

readImage :: BS.ByteString -> Int64 -> Vector R
readImage bytes n = vector $ fmap (fromIntegral . BS.index bytes . (n*28^2 + 16 +)) [0..783]

readLabel :: BS.ByteString -> Int64 -> Vector R
readLabel bytes n = vector [if x == i then 1 else 0 | x <- [0..9]] where
                               i = fromIntegral $ BS.index bytes (n + 8)

---------------------- Training ------------------------

train :: Network-> Int -> StdGen-> [(Vector R, Vector R)] -> Network
train net batchSize gen input = net' where
  net' = foldl' (foldl' processTrainingData) net batches
  batches = chunksOf batchSize $ shuffledInput
  shuffledInput = shuffle' input (length input) gen

processTrainingData :: Network -> (Vector R, Vector R) -> Network
processTrainingData net input = Network (layers net) bs ws (rate net) where
  ws = zipWith add (weights net) $ fmap (scale (-eta / l)) dw
  bs = zipWith add (biases net) $ fmap (scale (-eta / l)) db
  (dw, db) = unzip $ backprop net input
  l = fromIntegral $ length input
  eta = rate net

-- Returns nablaWeights and nableBiases in a tuple
backprop :: Network -> (Vector R, Vector R) -> [(Matrix R, Vector R)]
backprop net (y, x) = undefined where

  -- Propagate backwards calculating
  -- feltas thourhgout network, with input layer as head
  ds = foldl' calcDelta [outputDelta] list where
    outputDelta = cost' (head as) y * sigmoid' (head zs)
    list = init $ zip3 ws zs $ drop 1 as
    ws = reverse $ weights net
    calcDelta :: [Vector R] -> (Matrix R, Vector R, Vector R) -> [Vector R]
    calcDelta (d:ds) (w, z, a) = d':d:ds where
      d' = (tr' w #> d * sigmoid' z)

  -- Feed forward through the network
  -- calculating all zs and as, with the output layer as head
  (zs, as) = feedForward net x

  {-
  [(dw, db)] = foldl' calcDeltas [(calcDw d as, d)] $ init $ zip3 ws zs $ drop 1 as
ws = reverse $ weights net

calcDw :: Vector R -> [Vector R] -> R
  calcDw d (a:a':_) = d `dot` a'
  calcD (a:_) y = sigmoid' a * cost' a y

  calcDeltas :: [Vector R] -> (Matrix R, Vector R, Vector R) -> [Vector R]
  calcDeltas (d:ds) (w, z, a) = d':d:ds where
    d' = (w #> d * sigmoid' z)

  -- Calculate the deltas, with the input layer at index 0
  ds = foldl' calcDelta [(calcDw d as, delta as y)] $ init $ zip3 ws zs $ drop 1 as where
    calcDelta :: [Vector R] -> (Matrix R, Vector R, Vector R) -> [Vector R]
    calcDelta (d:ds) (w, z, a) = (w #> d * sigmoid' z):d:ds
    delta as y = sigmoid' a * cost' a y
    ws = reverse $ weights net
-}

-- Returns a tuple of zs and as throughout the network
-- with the last activation first (like a stack)
feedForward :: Network -> Vector R -> ([Vector R], [Vector R])
feedForward net a = foldl' ffd ([a], []) $ zip (weights net) (biases net)  where
  ffd :: ([Vector R], [Vector R]) -> (Matrix R, Vector R) -> ([Vector R], [Vector R])
  ffd (a:as, zs) (w, b) = (a':a:as, z:zs) where
    a' = sigmoid z
    z  = w #> a + b

{-
-- Perform stochastic gradient descent on the batch
-- and return a new network with updated weights and biases
processBatch :: Network -> [(Vector R, Vector R)] -> Network
processBatch net batch = foldl' processTrainingData net batch where
  processTrainingData :: Network -> (Vector R, Vector R) -> Network
  processTrainingData = \net input ->
    Network (layers net) bs ws (rate net) where
      ws = zipWith add (weights net) $ fmap (scale (-eta / l)) dw
      bs = zipWith add (biases net) $ fmap (scale (-eta / l)) db
      (dw, db) = unzip $ backprop net input
      l = fromIntegral $ length batch
      eta = rate net
-}

-- Returns a tuple of zs and as throughout the network
-- with the last activation first (like a stack)
{-
oldfeedForward :: Network -> Vector R -> ([Vector R], [Vector R])
oldfeedForward net a = ffd (weights net) (biases net) ([a], []) where
  ffd :: [Matrix R] -> [Vector R] -> ([Vector R], [Vector R]) -> ([Vector R], [Vector R])
  ffd [] [] a = a
  ffd (w:ws) (b:bs) (a:as, zs) = ffd ws bs (a':a:as, z:zs) where
    a' = sigmoid z
    z  = w #> a + b
-}
  -- Do the backpropagation
--forM ds $ \(lbl, img) ->
--                                    let a = feedForward n img

-- Cost function
cost w b n xs a = 1/(2*n) * (sum [(norm_2 (y_ x - a))^2 | x <- xs])

-- Derivative of cost
cost' :: Vector R -> Vector R -> Vector R
cost' a y = a - y

sigmoid :: Floating a => a -> a
sigmoid z = 1/(1 + exp(z))

sigmoid' :: Floating a => a -> a
sigmoid' z = sigmoid(z)*(1 - sigmoid(z))

--------------------  Utility   ---------------------------

randomMatrix :: Int -> Int -> Int -> Matrix R
randomMatrix x y seed = (x><y) . toList . randomVector seed Gaussian $ x * y

y_ :: Vector R -> Vector R
y_ = undefined

-- Test stuff
{-
getImage s n = fromIntegral . BS.index s . (n*28^2 + 16 +) <$> [0..28^2 - 1]

render n = let s = " .:oO@" in s !! (fromIntegral n * length s `div` 256)
pasta = do
  s <- decompress <$> BS.readFile "images.gz"
  l <- decompress <$> BS.readFile "labels.gz"
  let n = 1

  -- Labels 8 bytes
  -- Images has offset 16 8 * 28  28?

  --n <- (`mod` 60000) <$> randomIO
  let t =  (getImage s n :: [Double])
  putStrLn $ show t
  {-
  putStr . unlines $
    [(render . BS.index s . (n*28^2 + 16 + r*28 +)) <$> [0..27] | r <- [0..27]]
  print $ BS.index l (n + 8)
-}
-}

max_stars_repo_name: SebastianCallh/haskell-neural-net

max_issues_repo_name: SebastianCallh/haskell-neural-net

max_forks_repo_name: SebastianCallh/haskell-neural-net


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
#ifndef MIN_VERSION_base
#define MIN_VERSION_base(x,y,z) 1
#endif
-----------------------------------------------------------------------------
-- |
-- Module      :  Data.Complex.Lens
-- Copyright   :  (C) 2012-13 Edward Kmett
-- License     :  BSD-style (see the file LICENSE)
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  experimental
-- Portability :  non-portable
--
-- Lenses and traversals for complex numbers
--
----------------------------------------------------------------------------
module Data.Complex.Lens
  ( _realPart
  , _imagPart
  , _polar
  , _magnitude
  , _phase
  , _conjugate
  ) where

import Control.Applicative
import Control.Lens
import Data.Complex

-- $setup
-- >>> import Debug.SimpleReflect
-- >>> let { a ≈ b = abs (a - b) < 1e-6; infix 4 ≈ }

-- | Access the 'realPart' of a 'Complex' number.
--
-- >>> (a :+ b)^._realPart
-- a
--
-- >>> a :+ b & _realPart *~ 2
-- a * 2 :+ b
--
-- @'_realPart' :: 'Functor' f => (a -> f a) -> 'Complex' a -> f ('Complex' a)@
#if MIN_VERSION_base(4,4,0)
_realPart :: Lens' (Complex a) a
#else
_realPart :: RealFloat a => Lens' (Complex a) a
#endif
_realPart f (a :+ b) = (:+ b) <$> f a
{-# INLINE _realPart #-}

-- | Access the 'imagPart' of a 'Complex' number.
--
-- >>> (a :+ b)^._imagPart
-- b
--
-- >>> a :+ b & _imagPart *~ 2
-- a :+ b * 2
--
-- @'_imagPart' :: 'Functor' f => (a -> f a) -> 'Complex' a -> f ('Complex' a)@
#if MIN_VERSION_base(4,4,0)
_imagPart :: Lens' (Complex a) a
#else
_imagPart :: RealFloat a => Lens' (Complex a) a
#endif
_imagPart f (a :+ b) = (a :+) <$> f b
{-# INLINE _imagPart #-}

-- | This isn't /quite/ a legal 'Lens'. Notably the
--
-- @'view' l ('set' l b a) = b@
--
-- law is violated when you set a 'polar' value with 0 'magnitude' and non-zero
-- 'phase' as the 'phase' information is lost, or with a negative 'magnitude'
-- which flips the 'phase' and retains a positive 'magnitude'. So don't do
-- that!
--
-- Otherwise, this is a perfectly cromulent 'Lens'.
_polar :: RealFloat a => Iso' (Complex a) (a,a)
_polar = iso polar (uncurry mkPolar)
{-# INLINE _polar #-}

-- | Access the 'magnitude' of a 'Complex' number.
--
-- >>> (10.0 :+ 20.0) & _magnitude *~ 2
-- 20.0 :+ 40.0
--
-- This isn't /quite/ a legal 'Lens'. Notably the
--
-- @'view' l ('set' l b a) = b@
--
-- law is violated when you set a negative 'magnitude'. This flips the 'phase'
-- and retains a positive 'magnitude'. So don't do that!
--
-- Otherwise, this is a perfectly cromulent 'Lens'.
--
-- Setting the 'magnitude' of a zero 'Complex' number assumes the 'phase' is 0.
_magnitude :: RealFloat a => Lens' (Complex a) a
_magnitude f c = setMag <$> f r
  where setMag r' | r /= 0    = c * (r' / r :+ 0)
                  | otherwise = r' :+ 0
        r = magnitude c
{-# INLINE _magnitude #-}

-- | Access the 'phase' of a 'Complex' number.
--
-- >>> (mkPolar 10 (2-pi) & _phase +~ pi & view _phase) ≈ 2
-- True
--
-- This isn't /quite/ a legal 'Lens'. Notably the
--
-- @'view' l ('set' l b a) = b@
--
-- law is violated when you set a 'phase' outside the range @(-'pi', 'pi']@.
-- The phase is always in that range when queried. So don't do that!
--
-- Otherwise, this is a perfectly cromulent 'Lens'.
_phase :: RealFloat a => Lens' (Complex a) a
_phase f c = setPhase <$> f theta
  where setPhase theta' = c * cis (theta' - theta)
        theta = phase c
{-# INLINE _phase #-}

-- | Access the 'conjugate' of a 'Complex' number.
--
-- >>> (2.0 :+ 3.0) & _conjugate . _imagPart -~ 1
-- 2.0 :+ 4.0
--
-- >>> (mkPolar 10.0 2.0 ^. _conjugate . _phase) ≈ (-2.0)
-- True
_conjugate :: RealFloat a => Iso' (Complex a) (Complex a)
_conjugate = involuted conjugate
{-# INLINE _conjugate #-}

max_stars_repo_name: haasn/lens

max_issues_repo_name: haasn/lens

max_forks_repo_name: haasn/lens


text:
{-# LANGUAGE FlexibleContexts #-}
-- |
-- Module    : Statistics.Autocorrelation
-- Copyright : (c) 2009 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Functions for computing autocovariance and autocorrelation of a
-- sample.

module Statistics.Autocorrelation
    (
      autocovariance
    , autocorrelation
    ) where

import Prelude hiding (sum)
import Statistics.Function (square)
import Statistics.Sample (mean)
import Statistics.Sample.Internal (sum)
import qualified Data.Vector.Generic as G

-- | Compute the autocovariance of a sample, i.e. the covariance of
-- the sample against a shifted version of itself.
autocovariance :: (G.Vector v Double, G.Vector v Int) => v Double -> v Double
autocovariance a = G.map f . G.enumFromTo 0 $ l-2
  where
    f k = sum (G.zipWith (*) (G.take (l-k) c) (G.slice k (l-k) c))
          / fromIntegral l
    c   = G.map (subtract (mean a)) a
    l   = G.length a

-- | Compute the autocorrelation function of a sample, and the upper
-- and lower bounds of confidence intervals for each element.
--
-- /Note/: The calculation of the 95% confidence interval assumes a
-- stationary Gaussian process.
autocorrelation :: (G.Vector v Double, G.Vector v Int) => v Double -> (v Double, v Double, v Double)
autocorrelation a = (r, ci (-), ci (+))
  where
    r           = G.map (/ G.head c) c
      where c   = autocovariance a
    dllse       = G.map f . G.scanl1 (+) . G.map square $ r
      where f v = 1.96 * sqrt ((v * 2 + 1) / l)
    l           = fromIntegral (G.length a)
    ci f        = G.cons 1 . G.tail . G.map (f (-1/l)) $ dllse

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
-- |
-- Module      : Southpaw.Picasso.Shapes
-- Description : Various functions for creating mathematical shapes
-- Copyright   : (c) Jonatan H Sundqvist, year
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created date year

-- TODO | -
--        -

-- SPEC | -
--        -

---------------------------------------------------------------------------------------------------
-- API
---------------------------------------------------------------------------------------------------
module Southpaw.Picasso.Shapes where

---------------------------------------------------------------------------------------------------
-- We'll need these
---------------------------------------------------------------------------------------------------
import Data.Complex

import Southpaw.Math.Constants

---------------------------------------------------------------------------------------------------
-- Function
---------------------------------------------------------------------------------------------------
-- Geometry ---------------------------------------------------------------------------------------
-- |
-- TODO: Start angle
-- TODO: Invalid arguments (eg. sides < 3) (use Maybe?)
-- TODO: Make polymorphic
-- TODO: Move to geometry section
-- polygon :: (Floating f, RealFloat f) => Int -> f -> Complex f -> [Complex f]
polygon :: Integral int => int -> Double -> Complex Double -> [Complex Double]
polygon sides radius origin = [ let θ = arg n in origin + ((radius * cos θ):+(radius * sin θ)) | n <- [1..sides]]
    where arg n = fromIntegral n * (2*π)/fromIntegral sides

-- |
-- TODO: Simplify and expound comments
arrow :: Complex Double -> Complex Double -> Double -> Double -> Double -> [Complex Double]
arrow from to sl sw hw = [from     + straight (sw/2), --
                          shaftEnd + straight (sw/2), --
                          shaftEnd + straight (hw/2), --
                          to,                         --
                          shaftEnd - straight (hw/2), --
                          shaftEnd - straight (sw/2), --
                          from     - straight (sw/2)] --
    where along a b distance = (a +) . mkPolar distance . snd . polar $ b-a -- Walk distance along the from a to b
          normal a b = let (mag, arg) = polar (b-a) in mkPolar mag (arg+π/2)
          shaftEnd = along from to sl --
          straight = along (0:+0) (normal from to) -- Vector perpendicular to the centre line


text:
-- |
-- Module:       Math.NumberTheory.RootsOfUnityTests
-- Copyright:    (c) 2018 Bhavik Mehta
-- License:      MIT
-- Maintainer:   Andrew Lelechenko <andrew.lelechenko@gmail.com>
--
-- Tests for Math.NumberTheory.RootsOfUnity
--

module Math.NumberTheory.RootsOfUnityTests where

import Test.Tasty

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

import Math.NumberTheory.DirichletCharacters (toRootOfUnity, toComplex)
import Math.NumberTheory.TestUtils (testSmallAndQuick, Positive(..))

rootOfUnityTest :: Integer -> Positive Integer -> Bool
rootOfUnityTest n (Positive d) = toComplex ((d `div` gcd n d) `stimes` toRootOfUnity (n % d)) == (1 :: Complex Double)

testSuite :: TestTree
testSuite = testSmallAndQuick "RootOfUnity contains roots of unity" rootOfUnityTest

max_stars_repo_name: cjlarose/arithmoi

max_issues_repo_name: cjlarose/arithmoi

max_forks_repo_name: cjlarose/arithmoi


text:
{-# language TypeOperators #-}
{-# language FlexibleInstances #-}
{-# language UndecidableInstances #-}
{-# language MonoLocalBinds #-}
module Data.Functor.Comm
  ( Comm
  , CommBind
  ) where

import Control.Applicative (ZipList, Const)
import Control.Applicative.Backwards
import Control.Commutative
import Control.Comonad
import Control.Comonad.Env
import Control.Comonad.Traced
import Control.Monad.Trans.Identity
import Control.Monad.Trans.Maybe
import Control.Monad.Trans.Reader
import Control.Monad.Trans.Writer.Lazy as Lazy
import Control.Monad.Trans.Writer.Strict as Strict
import Data.Complex
import Data.Functor.Bind
import Data.Functor.Compose
import Data.Functor.Identity
import Data.Functor.Product as Functor
import Data.Functor.Reverse
import Data.Map.Lazy
import Data.Monoid as Monoid
import Data.Proxy
import Data.Semigroup as Semigroup
import Data.Semigroup.Commutative
import Data.Tagged
import GHC.Generics as Generics
import Linear
import Linear.Plucker
import Linear.V

-- |
-- @
-- ('<.>') = 'flip' ('<..>')
-- @
class Apply f => Comm f

instance (Comm f, Comm g) => Comm (Compose f g)
instance (Comm m, CommutativeSemigroup w) => Comm (Strict.WriterT w m)
instance (Comm m, CommutativeSemigroup w) => Comm (Lazy.WriterT w m)
instance CommutativeSemigroup a => Comm ((,) a)

-- instance (CommutativeSemigroup a, CommutativeSemigroup b) => Comm ((,,) a b)
-- instance (CommutativeSemigroup a, CommutativeSemigroup b, CommutativeSemigroup c) => Comm ((,,,) a b c)
-- instance (CommutativeSemigroup a, CommutativeSemigroup b, CommutativeSemigroup c, CommutativeSemigroup d) => Comm ((,,,,) a b c d)

instance Comm ((->) e)
instance Comm m => Comm (ReaderT e m)
instance Comm V4
instance Comm V3
instance Comm V2
instance Comm Linear.V1
instance Comm V0
instance Comm Plucker
instance Dim n => Comm (V n)
instance Comm Complex
instance Comm Quaternion
instance Comm Identity
instance Comm m => Comm (IdentityT m)
instance (Comm m, Comm n) => Comm (Functor.Product m n)
instance CommutativeMonad m => Comm (MaybeT m)
instance Comm m => Comm (Reverse m)
instance Comm m => Comm (Backwards m)
instance CommutativeSemigroup a => Comm (Const a)
instance Comm (Tagged a)
instance Comm Proxy
-- instance Comm m => Comm (Kleisli m a)
instance Comm (Cokleisli w a)
instance Comm ZipList
instance Comm Option -- NB: Not CommutativeSemigroup!
-- semigroupoids missing instances:
instance Comm f => Comm (Alt f)
instance Comm Semigroup.Last -- NB: Not CommutativeSemigroup!
instance Comm Semigroup.First -- NB: Not CommutativeSemigroup!
instance Comm Monoid.Last -- NB: Not CommutativeSemigroup!
instance Comm Monoid.First -- NB: Not CommutativeSemigroup!
instance Comm Monoid.Product -- NB: Not CommutativeSemigroup!
instance Comm Monoid.Sum -- NB: Not CommutativeSemigroup!

instance Comm m => Comm (M1 i c m)
instance (Comm f, Comm g) => Comm (f :.: g)
instance (Comm f, Comm g) => Comm (f :*: g)
instance CommutativeSemigroup c => Comm (K1 i c) -- missing in base
instance Comm f => Comm (Rec1 f)
instance Comm Par1
instance Comm U1
instance (Comm w, CommutativeSemigroup e) => Comm (EnvT e w)
instance Comm w => Comm (TracedT e w)
-- instance Comm f => Comm (Ap f) -- ghc 8.6

class    (Comm m, Bind m) => CommBind m
instance (Comm m, Bind m) => CommBind m

-- instances that don't extend to Commutative
instance Ord k => Comm (Map k)
instance Comm Generics.V1

max_stars_repo_name: ekmett/abelian

max_issues_repo_name: ekmett/abelian

max_forks_repo_name: ekmett/abelian


text:
{-# LANGUAGE RecordWildCards #-}

module HEP.Kinematics.Variable.M2
    (
      M2Solution (..)

    , mkInput

    , m2XXSQP
    , m2CXSQP
    , m2XCSQP
    , m2CCSQP

    , m2XXAugLag
    , m2CXAugLag
    , m2XCAugLag
    , m2CCAugLag
    ) where

import HEP.Kinematics
import HEP.Kinematics.Vector.LorentzVector (setXYZT)

import Numeric.LinearAlgebra               (Vector, fromList, toList)
import Numeric.NLOPT

-- import Numeric.GSL.Differentiation
-- import Debug.Trace

-- | All the components are rescaled by '_scale'.
data InputKinematics = InputKinematics
                       { _p1     :: FourMomentum  -- ^ p1 = a1 + b1
                       , _p2     :: FourMomentum  -- ^ p2 = a2 + b2
                       , _q1     :: FourMomentum  -- ^ q1 = b1
                       , _q2     :: FourMomentum  -- ^ q2 = b2
                       , _ptmiss :: TransverseMomentum
                         -- | squared mass of the invisible particle
                       , _mInvSq :: !Double
                       , _scale  :: !Double }

-- | creates the input for the process of
--
--   A1 + A2 --> a1 B1 + a2 B2 --> a1 b1 C1 + a2 b2 C2.
--
mkInput :: [FourMomentum]      -- ^ [a1, a2]
        -> [FourMomentum]      -- ^ [b1, b2]
        -> TransverseMomentum  -- ^ pT(miss)
        -> Double              -- ^ M_{invisible}
        -> Maybe InputKinematics
mkInput as bs ptmiss mInv
    | length as /= 2 || length bs /= 2 = Nothing
    | otherwise                        = do
          let ps = zipWith (+) as bs
              p1 = head ps
              p2 = ps !! 1

              m1 = mass p1
              m2 = mass p2
              mInvSq = mInv * mInv
              scaleSq = m1 * m1 + m2 * m2 + 2 * mInvSq
          if scaleSq <= 0  -- why?
              then Nothing
              else do let scale = sqrt scaleSq
                          p1'     = p1        ^/ scale
                          p2'     = p2        ^/ scale
                          q1'     = head bs   ^/ scale
                          q2'     = (bs !! 1) ^/ scale
                          ptmiss' = ptmiss    ^/ scale
                          mInv'   = mInv       / scale
                      return $ InputKinematics { _p1     = p1'
                                               , _p2     = p2'
                                               , _q1     = q1'
                                               , _q2     = q2'
                                               , _ptmiss = ptmiss'
                                               , _mInvSq = mInv' * mInv'
                                               , _scale  = scale }

data M2Solution = M2Solution { _M2    :: Double
                             , _k1sol :: FourMomentum
                             , _k2sol :: FourMomentum
                             } deriving Show

getM2Solution :: InputKinematics -> Either Result Solution -> Maybe M2Solution
getM2Solution inp@InputKinematics {..} sol =
    case sol of
        Left _                   -> Nothing
        -- sol0@(Right (Solution m2 ks _)) -> do
        Right (Solution m2 ks _) -> do
            let Invisibles k1 k2 = mkInvisibles inp (vecToVars ks)
            -- traceM ("sol = " ++ show sol0)
            return $ M2Solution { _M2    = m2  * _scale
                                , _k1sol = k1 ^* _scale
                                , _k2sol = k2 ^* _scale }

{-
initialGuess :: InputKinematics -> Vector Double
initialGuess InputKinematics { _ptmiss = ptmiss } =
    fromList [0.5 * px ptmiss, 0.5 * py ptmiss, 0, 0]
-}
initialGuess :: InputKinematics -> Vector Double
initialGuess inp@InputKinematics {..} =
    let objf ks = let Invisibles k1 k2 = mkInvisibles inp (vecToVars ks)
                  in invariantMass [_p1, _p2, k1, k2]
        ks0 = fromList [0.5 * px _ptmiss, 0.5 * py _ptmiss, 0, 0]
        eps' = eps * 0.1
        stop = MaximumEvaluations 1000
               :| [ ObjectiveAbsoluteTolerance eps'
                  , ObjectiveRelativeTolerance eps']
        algorithm = NELDERMEAD objf [] Nothing
        problem = LocalProblem 4 stop algorithm
        sol = minimizeLocal problem ks0
    in case sol of
        Left _                     -> ks0
        Right (Solution _ ksSol _) -> ksSol

eps :: Double
eps = 1e-3

stopObjCond :: [StoppingCondition]
stopObjCond = [ObjectiveAbsoluteTolerance eps', ObjectiveRelativeTolerance eps']
  where eps' = eps * 1e-3

type MultivarFunc = Vector Double -> (Double, Vector Double)

m2SQP :: [InputKinematics -> MultivarFunc]  -- ^ constraint functions
      -> Maybe InputKinematics
      -> Maybe M2Solution
m2SQP _   Nothing                         = Nothing
m2SQP cfs (Just inp@InputKinematics {..}) =
    let -- objective function with gradient
        objfD = m2ObjF inp

        -- constraint function with gradient
        constraints' = ($ inp) <$> cfs
        constraints  = (\cf -> EqualityConstraint (Scalar cf) eps)
                       <$> constraints'

        stop = MaximumEvaluations 1000 :| stopObjCond
        algorithm = SLSQP objfD [] [] constraints
        problem = LocalProblem 4 stop algorithm

        sol = minimizeLocal problem (initialGuess inp)
    in getM2Solution inp sol

m2XXSQP, m2CXSQP, m2XCSQP, m2CCSQP :: Maybe InputKinematics -> Maybe M2Solution
m2XXSQP = m2SQP []
m2CXSQP = m2SQP [constraintA]
m2XCSQP = m2SQP [constraintB]
m2CCSQP = m2SQP [constraintA, constraintB]

m2AugLag :: [InputKinematics -> MultivarFunc]  -- ^ constraint functions
         -> Maybe InputKinematics
         -> Maybe M2Solution
m2AugLag _   Nothing                         = Nothing
m2AugLag cfs (Just inp@InputKinematics {..}) =
    let -- objective function without gradient
        -- objf  = m2ObjF inp
        objf  = getResultOnly (m2ObjF inp)

        -- constraint function without gradient
        constraints' = ($ inp) <$> cfs
        -- constraints  = (\cf -> EqualityConstraint (Scalar cf) eps)
        --                <$> constraints'
        constraints  = (\cf -> EqualityConstraint (Scalar (getResultOnly cf))
                               eps) <$> constraints'

        stop = MaximumEvaluations 5000 :| stopObjCond
        -- algorithm = VAR1 objf Nothing
        -- algorithm = LBFGS objf Nothing
        -- algorithm = SLSQP objf [] [] constraints
        algorithm = NELDERMEAD objf [] Nothing

        subproblem = LocalProblem 4 stop algorithm
        -- problem = AugLagProblem [] constraints (AUGLAG_EQ_LOCAL subproblem)
        problem = AugLagProblem constraints [] (AUGLAG_EQ_LOCAL subproblem)

        sol = minimizeAugLag problem (initialGuess inp)
    in getM2Solution inp sol
  where
    getResultOnly :: MultivarFunc -> Vector Double -> Double
    getResultOnly fD ks = let (result, _) = fD ks in result

m2XXAugLag, m2CXAugLag, m2XCAugLag, m2CCAugLag
    :: Maybe InputKinematics -> Maybe M2Solution
m2XXAugLag = m2AugLag []
m2CXAugLag = m2AugLag [constraintA]
m2XCAugLag = m2AugLag [constraintB]
m2CCAugLag = m2AugLag [constraintA, constraintB]

-- | (k1x, k1y, k1z, k2z).
type Variables = (Double, Double, Double, Double)

-- | unsafe transformation, but it would be fine.
vecToVars :: Vector Double -> Variables
vecToVars ks = let [k1x, k1y, k1z, k2z] = toList ks in (k1x, k1y, k1z, k2z)

data Invisibles = Invisibles { _k1 :: FourMomentum , _k2 :: FourMomentum }

mkInvisibles :: InputKinematics -> Variables -> Invisibles
mkInvisibles InputKinematics {..} (k1x, k1y, k1z, k2z) = Invisibles k1 k2
  where
    eInv1 = sqrt (k1x * k1x + k1y * k1y + k1z * k1z + _mInvSq)
    k1 = eInv1 `seq` setXYZT k1x k1y k1z eInv1

    k2x = px _ptmiss - k1x
    k2y = py _ptmiss - k1y
    eInv2 = sqrt (k2x * k2x + k2y * k2y + k2z * k2z + _mInvSq)
    k2 = eInv2 `seq` setXYZT k2x k2y k2z eInv2

-- | the unknowns are (k1x, k2x, k1z, k2z).
m2ObjF :: InputKinematics -> MultivarFunc
m2ObjF inp@InputKinematics {..} ks =
    if m1 < m2 then (m2, grad2) else (m1, grad1)
  where
    invs = mkInvisibles inp (vecToVars ks)
    (grad1, grad2, m1, m2) = m2Grad inp invs _p1 _p2 ks

m2Grad :: InputKinematics
       -> Invisibles
       -> FourMomentum  -- ^ p1 (or q1)
       -> FourMomentum  -- ^ p2 (or q2)
       -> Vector Double
       -> (Vector Double, Vector Double, Double, Double)
m2Grad InputKinematics {..} (Invisibles k1 k2) p1 p2 ks = (d1, d2, m1, m2)
  where
    (k1x, k1y, k1z, k2z) = vecToVars ks

    (e1, p1x, p1y, p1z) = epxpypz p1
    m1  = invariantMass [p1, k1]
    m1' = safeDivisor m1
    r1 = e1 / safeDivisor (energy k1)
    d1 = fromList $ ( / m1') <$> [ r1 * k1x - p1x
                                 , r1 * k1y - p1y
                                 , r1 * k1z - p1z
                                 , 0 ]

    (e2, p2x, p2y, p2z) = epxpypz p2
    m2  = invariantMass [p2, k2]
    m2' = safeDivisor m2
    r2 = e2 / safeDivisor (energy k2)
    d2 = fromList $ ( / m2') <$> [ r2 * (k1x - px _ptmiss) + p2x
                                 , r2 * (k1y - py _ptmiss) + p2y
                                 , 0
                                 , r2 * k2z - p2z ]

constraintF :: FourMomentum -> FourMomentum -> InputKinematics -> MultivarFunc
constraintF p1 p2 inp ks = (m1 - m2, grad1 - grad2)
  where
    invs = mkInvisibles inp (vecToVars ks)
    (grad1, grad2, m1, m2) = m2Grad inp invs p1 p2 ks

constraintA, constraintB :: InputKinematics -> MultivarFunc
constraintA inp@InputKinematics {..} = constraintF _p1 _p2 inp
constraintB inp@InputKinematics {..} = constraintF _q1 _q2 inp

safeDivisor :: Double -> Double
safeDivisor x | x >= 0    = max eps0 x
              | otherwise = min (-eps0) x
  where eps0 = 1.0e-8

{-
m2ObjF :: InputKinematics -> MultivarFunc
m2ObjF inp ks = (m, d)
  where
    ks0 = vecToVars ks
    m = fst $ m2ObjF' inp ks0
    d = m2ObjFDNumeric inp ks0

m2ObjF' :: InputKinematics -> Variables -> (Double, Vector Double)
m2ObjF' inp@InputKinematics {..} ks =
    if m1 < m2 then (m2, grad2) else (m1, grad1)
  where
    invs@(Invisibles k1 k2) = mkInvisibles inp ks
    m1 = invariantMass [_p1, k1]
    m2 = invariantMass [_p2, k2]
    (k1x, k1y, k1z, k2z) = ks
    (grad1, grad2, _) = m2Grad inp invs _p1 _p2 (fromList [k1x, k1y, k1z, k2z])

m2ObjFDNumeric :: InputKinematics -> Variables -> Vector Double
m2ObjFDNumeric inp@InputKinematics {..} = fDNumeric f
  where f ks = fst (m2ObjF' inp ks)

m2GradNumeric :: InputKinematics
              -> Invisibles
              -> FourMomentum  -- ^ p1 (or q1)
              -> FourMomentum  -- ^ p2 (or q2)
              -> Vector Double
              -> (Vector Double, Vector Double, Double)
m2GradNumeric inp@(InputKinematics {..}) _ p1 p2 ks = (d1, d2, m1 - m2)
  where
    ks0 = vecToVars ks

    m1F ks' = let Invisibles k1 _ = mkInvisibles inp ks'
              in invariantMass [p1, k1]
    m1 = m1F ks0
    d1 = fDNumeric m1F ks0

    m2F ks' = let Invisibles _ k2 = mkInvisibles inp ks'
              in invariantMass [p2, k2]
    m2 = m2F ks0
    d2 = fDNumeric m2F ks0

fDNumeric :: (Variables -> Double) -> Variables -> Vector Double
fDNumeric f (k1x, k1y, k1z, k2z) = fromList [fk1xD, fk1yD, fk1zD, fk2zD]
  where
    fk1x a = f (a, k1y, k1z, k2z)
    fk1y a = f (k1x, a, k1z, k2z)
    fk1z a = f (k1x, k1y, a, k2z)
    fk2z a = f (k1x, k1y, k1z, a)

    epsN = 1e-8

    fk1xD = fst $ derivCentral epsN fk1x k1x
    fk1yD = fst $ derivCentral epsN fk1y k1y
    fk1zD = fst $ derivCentral epsN fk1z k1z
    fk2zD = fst $ derivCentral epsN fk2z k2z
-}


text:
{-|
Module      : Distributions
Description : Definition of the Distributions used in Models and associated functions.
Copyright   : (c) Julian Kopka Larsen, 2015
Stability   : experimental
-}
module Distributions (
    lnormalInvWishart,
    lnormalInvWishartSS,
    lMixNormWish,
    delta,
    Expr,
    dlNormW,
    lNormW
) where

import Math (X,
             mlgamma,
             scatterMatrix,
             meanv,
             lnDeterminant)

import Numeric.LinearAlgebra (Matrix,
                              Vector,
                              (<>),
                              ident,
                              dim,
                              scale,
                              trans,
                              takeDiag,
                              outer)

import Numeric.LinearAlgebra.Data (konst, size)

import Partition (Partition, Component, groups, group)

import LikSpecLang

-- | Delta log Normal-Wishard likelihood for a mixture of components.
dlNormW :: X -> (Partition -> Partition -> Double)
dlNormW x = delta (lMixNormWish x)

-- | Log Normal-Wishard likelihood for a mixture of components.
lNormW :: X -> (Partition -> Double)
lNormW x = compile (lMixNormWish x)

-- | Log Normal-Wishard likelihood for a mixture of components. Same as 'lNormW' but defined without the "LikSpecLang" (LSL)
lMixNormWish ::X -> Expr Partition Double
lMixNormWish x = Mixture (Fun lnormalInvWishart) % selection x -- The '%' operator "pipes" expressions together

selection :: X -> Expr Partition [X]
selection x = Fun (groups x)

-- | Log Normal-Wishard likelihood for a single component
lnormalInvWishartSS :: Int -> Matrix Double -> Vector Double-> Double
lnormalInvWishartSS n s mu = - fromIntegral n * d * 2 * log pi
                      + mlgamma d (vn/2)
                      - mlgamma d (v0/2)
                      + (v0/2) * ldalpha0
                      - (vn/2) * ldalphaN
                      + (d/2) * ((log k0)-(log kn))
                    where   d = fromIntegral $ fst $ size s :: Double
                            v0 = 100 -- hyperparameter must be (v0 > d)
                            vn = v0 + fromIntegral n
                            alpha0 = ident $ fst $ size s -- hyperparameter for the shape of the covariance matrix
                            alphaN = alpha0 + s + scale ((k0*fromIntegral n)/(k0+kn)) ((mu-mu0) `outer` (mu-mu0))
                            ldalpha0 = lnDeterminant alpha0
                            ldalphaN = lnDeterminant alphaN
                            k0 = 0.2 -- hyperparameter
                            kn = k0 + fromIntegral n
                            mu0 = (konst 0.0 $ fst $ size s ):: Vector Double -- hyperparameter

-- | Log Normal-Wishard likelihood for a single component
lnormalInvWishart :: X -> Double
lnormalInvWishart x = - n * d * 2 * log pi
                      + mlgamma d (vn/2)
                      - mlgamma d (v0/2)
                      + (v0/2) * ldalpha0
                      - (vn/2) * ldalphaN
                      + (d/2) * ((log k0)-(log kn))
                    where   n = fromIntegral $ length x
                            d = fromIntegral $ dim $ head x
                            v0 = 100 -- hyperparameter must be (v0 > d)
                            vn = v0 + n
                            alpha0 = ident $ dim $ head x -- hyperparameter for the shape of the covariance matrix
                            alphaN = alpha0 + s + scale ((k0*n)/(k0+kn)) ((xm-mu0) `outer` (xm-mu0))
                            ldalpha0 = lnDeterminant alpha0
                            ldalphaN = lnDeterminant alphaN
                            k0 = 0.2 -- hyperparameter
                            kn = k0 + n
                            s = scatterMatrix x
                            mu0 = (konst 0.0 $ dim $ head x ):: Vector Double -- hyperparameter
                            xm = meanv x

max_stars_repo_name: juliankopkalarsen/FpStats

max_issues_repo_name: juliankopkalarsen/FpStats

max_forks_repo_name: juliankopkalarsen/FpStats


text:
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE OverloadedStrings         #-}
module ML.Gym.Range
  ( GymRange(..)
  , combineRanges
  , getGymRange
  , gymRangeToDoubleLists
  , gymRangeToFloatLists
  ) where

import           ML.Gym.DType
import           ML.Gym.Util
import           ML.Gym.Value

import           Control.Applicative
import           Control.Arrow
import qualified Control.Exception        as E
import           Control.Monad            (join, void, (>=>))
import qualified CPython                  as Py
import qualified CPython.Constants        as Py
import qualified CPython.Protocols.Object as Py
import qualified CPython.Types            as Py
import qualified CPython.Types.Exception  as Py
import qualified CPython.Types.Module     as Py
import qualified Data.ByteString          as BS
import           Data.Complex             (Complex)
import           Data.Foldable            (toList)
import           Data.Maybe               (fromMaybe)
import qualified Data.Text                as T
import           System.Environment       (getArgs, getProgName)
import           System.Exit
import           System.IO
import           System.Random            (randomRIO)

import           Debug.Trace

data GymRange =
  GymRange GymDType
           (GymValue Double)    -- ^ Minimum values.
           (GymValue Double)    -- ^ Maximum values.
  deriving (Show, Eq)

gymRangeToDoubleLists :: GymRange -> ([Double],[Double])
gymRangeToDoubleLists (GymRange _ los his) = (gymValueTo1D los, gymValueTo1D his)

gymRangeToFloatLists :: GymRange -> ([Float],[Float])
gymRangeToFloatLists = (map realToFrac *** map realToFrac) . gymRangeToDoubleLists

gymRangeTypePrecedence :: GymRange -> Int
gymRangeTypePrecedence (GymRange dt _ _) = case dt of
  GymBool    -> 1
  GymInteger -> 2
  GymFloat   -> 3
  _          -> error $ "range type " ++ show dt ++ " not allowed"

gymRangeTypeFromPrecedence :: Int -> GymDType
gymRangeTypeFromPrecedence 1 = GymBool
gymRangeTypeFromPrecedence 2 = GymInteger
gymRangeTypeFromPrecedence 3 = GymFloat
gymRangeTypeFromPrecedence nr = error $ "gym Range with precedence " ++ show nr ++ " is unknown"

combineRanges :: [GymRange] -> GymRange
combineRanges [] = error "empty list ranges in maxRange"
-- combineRanges xs = uncurry (GymRange (gymRangeTypeFromPrecedence $ maximum (map gymRangeTypePrecedence xs))) (foldl1 mkLoHiRanges (map gymRangeToDoubleLists xs))
--   where mkLoHiRanges (los,his) (los', his') | length los' > length los = mkLoHiRanges (los ++ drop (length los) los', his ++ drop (length his) his') (los', his')
--                                             | length los > length los' = mkLoHiRanges (los,his) (los' ++ drop (length los') los, his' ++ drop (length his') his)
--                                             | otherwise = (zipWith min los los', zipWith max his his')

getGymRange :: Py.Object space => GymDType -> space -> IO (Maybe GymRange)
getGymRange tp@GymBool    = getGymRange' pyToBool fromBool (GymRange tp)
getGymRange tp@GymInteger = getGymRange' pyToInteger fromIntegral (GymRange tp)
getGymRange tp@GymFloat   = getGymRange' pyToDouble id (GymRange tp)
getGymRange tp@GymObject  = getGymRange' pyToDouble id (GymRange tp)
getGymRange tp            = error $ "cannot make range of type " ++ show tp

numPyArray :: Py.SomeObject -> IO Py.SomeObject
numPyArray obj = python $ Py.callMethodArgs obj "tolist" []

getGymRange' :: Py.Object space => (Py.SomeObject -> IO (Maybe a)) -> (a -> Double) -> (GymValue Double -> GymValue Double -> GymRange) -> space -> IO (Maybe GymRange)
getGymRange' convert mkGymVal gymRange space = do
  simple <- getLoHg (convert >=> return . fmap (GymScalar . mkGymVal)) gymRange space
  list <- getLoHg (\x -> fmap Gym1D <$> (numPyArray x >>= pyToListOf (fmap (fmap mkGymVal) . convert))) gymRange space
  ll <- getLoHg (\x -> fmap Gym2D <$> (numPyArray x >>= pyToListOfListOf (fmap (fmap mkGymVal) . convert))) gymRange space
  lll <- getLoHg (\x -> fmap Gym3D <$> (numPyArray x >>= pyToListOfListOfListOf (fmap (fmap mkGymVal) . convert))) gymRange space
  return $ simple <|> list <|> ll <|> lll
  where
    getLoHg convert gymRange space = do
      lo <- lowerBound space >>= convert
      hg <- upperBound space >>= convert
      return $ pure gymRange <*> lo <*> hg

mkSpaceRepr :: Py.Object self => self -> IO T.Text
mkSpaceRepr space = fromMaybe (error "could not get gym action space") <$> python (Py.callMethodArgs space "__repr__" [] >>= pyToText)

lowerBound :: Py.Object space => space -> IO Py.SomeObject
lowerBound space = python (Py.toUnicode "low" >>= Py.getAttribute space)

upperBound :: Py.Object space => space -> IO Py.SomeObject
upperBound space = python (Py.toUnicode "high" >>= Py.getAttribute space)

max_stars_repo_name: schnecki/gym-haskell

max_issues_repo_name: schnecki/gym-haskell

max_forks_repo_name: schnecki/gym-haskell


text:
{-# LANGUAGE OverloadedStrings #-}
module ML.Gym.Data
  ( GymData(..)
  , getGymData
  , gymObservationToDoubleList
  , gymObservationToFloatList
  ) where

import           ML.Gym.DType
import           ML.Gym.Range
import           ML.Gym.Space
import           ML.Gym.Util
import           ML.Gym.Value

import           Control.Applicative
import qualified Control.Exception        as E
import           Control.Monad            (void, zipWithM, (>=>))
import qualified CPython                  as Py
import qualified CPython.Constants        as Py
import qualified CPython.Protocols.Object as Py
import qualified CPython.Types            as Py
import qualified CPython.Types.Exception  as Py
import qualified CPython.Types.Module     as Py
import qualified Data.ByteString          as BS
import           Data.Complex             (Complex)
import           Data.Maybe               (fromMaybe)
import qualified Data.Text                as T
import           System.Environment       (getArgs, getProgName)
import           System.Exit
import           System.IO
import           System.Random            (randomRIO)

type GymData = GymValue Double

gymObservationToDoubleList :: GymData -> [Double]
gymObservationToDoubleList = gymValueTo1D

gymObservationToFloatList :: GymData -> [Float]
gymObservationToFloatList = map realToFrac . gymObservationToDoubleList

getGymData :: GymSpace -> Py.SomeObject -> IO (Maybe GymData)
getGymData (GymSpaceDiscrete nr) obj  = fmap (GymScalar . fromIntegral) <$> pyToInteger obj
getGymData (GymSpaceBox range@(GymRange dt _ _) _) obj  = getGymValueWithRange range obj
getGymData (GymSpaceTuple spaces) obj = do
  list <- fromMaybe (error "Gym discrete space 'n' not recognized!") <$> python (Py.toUnicode "spaces" >>= Py.getAttribute obj >>= pyToList)
  fmap GymTuple . sequence <$> zipWithM getGymData spaces list

fromScalar :: GymDType -> Py.SomeObject -> IO (Maybe (GymValue Double))
fromScalar GymBool obj    = fmap (GymScalar . fromBool) <$> pyToBool obj
fromScalar GymInteger obj = fmap (GymScalar . fromIntegral) <$> pyToInteger obj
fromScalar GymFloat obj   = fmap GymScalar <$> pyToDouble obj
fromScalar dt _           = error $ "missing implemenations in fromScalar: " ++ show dt

from1D :: GymDType -> Py.SomeObject -> IO (Maybe (GymValue Double))
from1D GymBool obj    = fmap (Gym1D . map fromBool) <$> pyToListOf pyToBool obj
from1D GymInteger obj = fmap (Gym1D . map fromIntegral) <$> pyToListOf pyToInteger obj
from1D GymFloat obj   = fmap Gym1D <$> pyToListOf pyToDouble obj
from1D dt _           = error $ "missing implemenations in fromScalar: " ++ show dt

from2D :: GymDType -> Py.SomeObject -> IO (Maybe (GymValue Double))
from2D GymBool obj    = fmap (Gym2D . map (map fromBool)) <$> pyToListOfListOf pyToBool obj
from2D GymInteger obj = fmap (Gym2D . map (map fromIntegral)) <$> pyToListOfListOf  pyToInteger obj
from2D GymFloat obj   = fmap Gym2D <$> pyToListOfListOf pyToDouble obj
from2D dt _           = error $ "missing implemenations in fromScalar: " ++ show dt

from3D :: GymDType -> Py.SomeObject -> IO (Maybe (GymValue Double))
from3D GymBool obj    = fmap (Gym3D . map (map (map fromBool))) <$> pyToListOfListOfListOf pyToBool obj
from3D GymInteger obj = fmap (Gym3D . map (map (map fromIntegral))) <$> pyToListOfListOfListOf  pyToInteger obj
from3D GymFloat obj   = fmap Gym3D <$> pyToListOfListOfListOf pyToDouble obj
from3D dt _           = error ("missing implementations in fromScalar: " ++ show dt)

getGymValueWithRange :: GymRange -> Py.SomeObject -> IO (Maybe (GymValue Double))
getGymValueWithRange (GymRange dt (GymScalar _) _) obj = fromScalar dt obj
getGymValueWithRange (GymRange dt (Gym1D _) _) obj     = numPyArray obj >>= from1D dt
getGymValueWithRange (GymRange dt (Gym2D _) _) obj     = numPyArray obj >>= from2D dt
getGymValueWithRange (GymRange dt (Gym3D _) _) obj     = numPyArray obj >>= from3D dt
getGymValueWithRange (GymRange dt (GymTuple xs) _) obj = do
  list <- fromMaybe (error "Gym discrete space 'n' not recognized!") <$> python (Py.toUnicode "spaces" >>= Py.getAttribute obj >>= pyToList)
  fmap GymTuple . sequence <$> zipWithM (\v o -> getGymValueWithRange (GymRange dt v (error "unused")) o) xs list

numPyArray :: Py.SomeObject -> IO Py.SomeObject
numPyArray obj = python $ Py.callMethodArgs obj "tolist" []


text:
module NeuralNetwork where

import           Numeric.LinearAlgebra as LA
import           Text.Printf ( printf )
import           Numeric.Morpheus.Activation ( relu, reluGradient, sigmoid,
                    sigmoidGradient, tanh_, tanhGradient)

data Activation = Relu | Sigmoid | Tanh | Id

-- Neural network layer: weights, biases, and activation
data Layer a = Layer (Matrix a) (Matrix a) Activation

type NeuralNetwork a = [Layer a]

-- | Weight and bias gradients
data Gradients a = Gradients (Matrix a) (Matrix a)

getActivation :: Activation -> (Matrix Double -> Matrix Double)
getActivation Id = id
getActivation Sigmoid = sigmoid
getActivation Relu = relu
getActivation Tanh = tanh_

getActivation' :: Activation -> Matrix R -> Matrix R -> Matrix R
getActivation' Id = flip const
getActivation' Sigmoid = \x dY -> sigmoidGradient x * dY
getActivation' Relu = \x dY -> reluGradient x * dY
getActivation' Tanh = \x dY -> tanhGradient x * dY

bias' :: Matrix R -> Matrix R
bias' dY = cmap (/ m) r
  where
    r = matrix (cols dY) $ map sumElements (toColumns dY)
    m = fromIntegral $ rows dY

linearW' :: (Numeric t, Fractional t) => Matrix t -> Matrix t -> Matrix t
linearW' x dy = cmap (/ m) (tr' x LA.<> dy)
  where
    m = fromIntegral $ rows x

linearX' :: Numeric t => Matrix t -> Matrix t -> Matrix t
linearX' w dy = dy LA.<> tr' w

-- | Gradient descent optimization
optimize :: R -> Int -> [Layer R] -> (Matrix R, Matrix R) -> [Layer R]
optimize lr iterN net0 dataSet = last $ take iterN (iterate step net0)
  where
    step net = zipWith f net dW
        where
            (_, dW) = pass net dataSet

    f :: Layer R -> Gradients R -> Layer R
    f (Layer w b act) (Gradients dW dB) =
        Layer (w - lr `scale` dW) (b - lr `scale` dB) act

pass :: NeuralNetwork R -> (Matrix R, Matrix R) -> (Matrix R, [Gradients R])
pass net (x, tgt) = (pred, grads)
  where
    (_, pred, grads) = _pass x net

    _pass inp [] = (loss', pred, [])
      where
        pred = sigmoid inp
        loss' = pred - tgt

    _pass inp (Layer w b sact:layers) = (dX, pred, Gradients dW dB:t)
      where
        lin = (inp LA.<> w) + b
        y = getActivation sact lin

        (dZ, pred, t) = _pass y layers

        dY = getActivation' sact lin dZ
        dW = linearW' inp dY
        dB = bias' dY
        dX = linearX' w dY

data AdamParameters = AdamParameters { _beta1 :: R, _beta2 :: R, _epsilon :: R, _lr :: R }

adamParams :: AdamParameters
adamParams = AdamParameters { _beta1 = 0.9
                            , _beta2 = 0.999
                            , _epsilon = 1e-8
                            , _lr = 0.001
                            }

-- | Adam optimization
optimizeAdam :: AdamParameters -> Int -> [Layer R] -> (Matrix R, Matrix R) -> [Layer R]
optimizeAdam p iterN w0 dataSet = w
  where
    s0 = map zf w0
    v0 = map zf w0
    zf (Layer a b _) = (zerosLike a, zerosLike b)
    zerosLike m = matrix c (replicate (r*c) 0.0)
      where
        r = rows m
        c = cols m
    (w, _, _) = adam p iterN (w0, s0, v0) dataSet

adam :: AdamParameters -> Int -> ([Layer R], [(Matrix R, Matrix R)], [(Matrix R, Matrix R)]) -> (Matrix R, Matrix R) -> ([Layer R], [(Matrix R, Matrix R)], [(Matrix R, Matrix R)])
adam p@AdamParameters { _lr = lr, _beta1 = beta1, _beta2 = beta2, _epsilon = epsilon } iterN (w0, s0, v0) dataSet = last $ take iterN (iterate step (w0, s0, v0))
      where
        step (w, s, v) = (wN, sN, vN)
          where
            (_, dW) = pass w dataSet

            sN = zipWith f2 s dW
            vN = zipWith f3 v dW
            wN = zipWith3 f w vN sN

        f (Layer w_ b_ sf) (vW, vB) (sW, sB) =
           Layer (w_ - lr `scale` vW / (sqrt sW `addC` epsilon))
                 (b_ - lr `scale` vB / (sqrt sB `addC` epsilon))
                 sf

        addC m c = cmap (+ c) m

        f2 (sW, sB) (Gradients dW dB) =
          ( beta2 `scale` sW + (1 - beta2) `scale` (dW^2)
          , beta2 `scale` sB + (1 - beta2) `scale` (dB^2))

        f3 (vW, vB) (Gradients dW dB) =
          ( beta1 `scale` vW + (1 - beta1) `scale` dW
          , beta1 `scale` vB + (1 - beta1) `scale` dB)

genWeights :: (Int, Int) -> IO (Matrix Double, Matrix Double)
genWeights (nin, nout) = do
  w <- _genWeights (nin, nout)
  b <- _genWeights (1, nout)
  return (w, b)
    where
      _genWeights (nin, nout) = do
          let k = sqrt (1.0 / fromIntegral nin)
          w <- randn nin nout
          return (k `scale` w)

makeCircles :: Int -> R -> R -> IO (Matrix R, Matrix R)
makeCircles m factor noise = do
  let rand' n = scale (2 * pi) <$> rand n 1
      m1 = m `div` 2
      m2 = m - (m `div` 2)

  r1 <- rand' m1
  r2 <- rand' m2
  ns <- scale noise <$> randn m 2

  let outerX = cos r1
      outerY = sin r1
      innerX = scale factor $ cos r2
      innerY = scale factor $ sin r2
      -- Merge them all
      x = fromBlocks [[outerX, outerY], [innerX, innerY]]

      -- Labels
      y1 = m1 >< 1 $ repeat 0
      y2 = m2 >< 1 $ repeat 1
      y = y1 === y2

  return (x + ns, y)

makeSpirals :: Int -> R -> IO (Matrix R, Matrix R)
makeSpirals m noise = do
  r0 <- scale (780 * 2 * pi / 360) . sqrt <$> rand m 1
  d1x0 <- scale noise <$> rand m 1
  d1y0 <- scale noise <$> rand m 1

  let d1x = d1x0 - cos r0 * r0
  let d1y = d1y0 + sin r0 * r0

  let x = fromBlocks [[d1x, d1y], [- d1x, - d1y]] / 10.0
  let y1 = m >< 1 $ repeat 0
  let y2 = m >< 1 $ repeat 1
  let y = y1 === y2
  return (x, y)

accuracy :: [Layer R] -> (Matrix R, Matrix R) -> R
accuracy net (dta, tgt) = 100 * (1 - e / m)
  where
    pred = net `inferBinary` dta
    e = sumElements $ abs (tgt - pred)
    m = fromIntegral $ rows tgt

inferBinary :: NeuralNetwork R -> Matrix R -> Matrix R
inferBinary net dta =
  let pred = forward net dta
  in cmap (\a -> if a < 0.5 then 0 else 1) pred

forward :: NeuralNetwork R -> Matrix R -> Matrix R
forward net dta = fst $ pass net (dta, undefined)

genNetwork :: [Int] -> [Activation] -> IO [Layer R]
genNetwork nodes activations = do
    weights <- mapM genWeights nodes'
    return (zipWith (\(w, b) a -> Layer w b a) weights activations)
  where
    nodes' = zip nodes (tail nodes)

experiment1 :: IO ()
experiment1 = do
  trainSet <- makeCircles 200 0.6 0.1
  testSet <- makeCircles 100 0.6 0.1

  let (dta, tgt) = trainSet

  (w1_rand, b1_rand) <- genWeights (2, 128)
  (w2_rand, b2_rand) <- genWeights (128, 1)

  let net = [ Layer w1_rand b1_rand Relu
            , Layer w2_rand b2_rand Id ]

  -- -- Alternatively:
  -- net0 <- genNetwork [2, 128, 1] [Relu, Id]

  let epochs = 1000
      lr = 0.001  -- Learning rate
      net' = optimize lr epochs net trainSet
      netA = optimizeAdam adamParams epochs net trainSet

  putStrLn $ printf "Circles problem, 1 hidden layer of 128 neurons, %d epochs" epochs
  putStrLn "---"

  putStrLn $ printf "Training accuracy (gradient descent) %.1f" (net' `accuracy` trainSet)
  putStrLn $ printf "Validation accuracy (gradient descent) %.1f\n" (net' `accuracy` testSet)

  putStrLn $ printf "Training accuracy (Adam) %.1f" (netA `accuracy` trainSet)
  putStrLn $ printf "Validation accuracy (Adam) %.1f\n" (netA `accuracy` testSet)

  putStrLn ""

experiment2 :: IO ()
experiment2 = do
  trainSet <- makeSpirals 200 0.5
  testSet <- makeSpirals 100 0.5

  let (dta, tgt) = trainSet
      epochs = 700

  putStrLn $ printf "Spirals problem, Adam, %d epochs" epochs
  putStrLn "---"
  putStrLn "1 hidden layer, 128 neurons (513 parameters)"
  net0 <- genNetwork [2, 128, 1] [Relu, Id]
  let net0' = optimizeAdam adamParams epochs net0 trainSet

  putStrLn $ printf "Training accuracy %.1f" (net0' `accuracy` trainSet)
  putStrLn $ printf "Validation accuracy %.1f\n" (net0' `accuracy` testSet)

  putStrLn "1 hidden layer, 512 neurons (2049 parameters)"
  net1 <- genNetwork [2, 512, 1] [Relu, Id]
  let net1' = optimizeAdam adamParams epochs net1 trainSet

  putStrLn $ printf "Training accuracy %.1f" (net1' `accuracy` trainSet)
  putStrLn $ printf "Validation accuracy %.1f\n" (net1' `accuracy` testSet)

  putStrLn "3 hidden layers, 40, 25, and 10 neurons (1416 parameters)"
  net2 <- genNetwork [2, 40, 25, 10, 1] [Relu, Relu, Relu, Id]
  let net2' = optimizeAdam adamParams epochs net2 trainSet

  putStrLn $ printf "Training accuracy %.1f" (net2' `accuracy` trainSet)
  putStrLn $ printf "Validation accuracy %.1f\n" (net2' `accuracy` testSet)

main :: IO ()
main = experiment1 >> experiment2


text:
module Regression where

import Control.Monad.State
import Numeric.LinearAlgebra                         (col, disp, Matrix, R, rand, scale, toLists, tr, Vector, (><))
import Numeric.LinearAlgebra.HMatrix                 (mul, sumElements)

meanMatrix :: Matrix R -> Int -> Double
meanMatrix inputs batchSize = (sumElements inputs) / (fromIntegral batchSize)

variance :: [Double] -> Double -> Double
variance inputs mean = sum [(x - mean)**2 | x <- inputs]

sigmoid :: Matrix R -> Matrix R
sigmoid inputs = 1.0 / (1.0 + exp (-1 * inputs))

forward :: Matrix R -> Matrix R -> Matrix R
forward inputs weights = sigmoid (mul inputs weights)

criterion :: Matrix R -> Matrix R -> Matrix R
criterion targets logits = (-1 * targets) * (log logits) - (1 - targets) * (log $ 1 - logits)

loss :: Matrix R -> Matrix R -> Int -> Double
loss targets logits batchSize = meanMatrix (criterion targets logits) batchSize

getGradient :: Matrix R -> Matrix R -> Matrix R -> Int -> Matrix R
getGradient inputs targets logits batchSize = (mul (tr inputs) (logits - targets)) / (fromIntegral batchSize)

updateWeights :: Double -> Matrix R -> Matrix R -> Matrix R
updateWeights learningRate gradient weights = (weights - (learningRate `scale` gradient))

getNewWeights :: Matrix R -> Matrix R -> Matrix R -> Int -> Matrix R
getNewWeights inputs targets weights batchSize =
    updateWeights 0.001 (getGradient inputs targets (forward inputs weights) batchSize) weights

runTrain :: Matrix R -> Matrix R -> Int -> State (Matrix R) (Matrix R)
runTrain inputs targets batchSize = do
    weights <- get
    forM_ [0..500] $ \e -> do
        put (getNewWeights inputs targets weights batchSize)
    return weights

activate :: [Double] -> [Double]
activate predictions = [if x > 0.5 then 1.0 else 0.0 | x <- predictions ]

max_stars_repo_name: staturecrane/haskell-ml

max_issues_repo_name: staturecrane/haskell-ml

max_forks_repo_name: staturecrane/haskell-ml


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}

module Elescore.Api.Data
  ( DataApi
  , dataServer
  ) where

import           ClassyPrelude                         hiding (Handler)
import           Data.DateTime                         (DateTime)
import qualified Data.Map                              as M
import qualified Data.Text                             as T
import           Servant
import           Servant.Pagination

import           Elescore.Api.Types
import           Elescore.IdTypes
import           Elescore.Projection.DisruptionsPerDay (DisruptionsPerDay)
import           Statistics.IQR                        (average)

type AppliedRanges ranges a = Ranges ranges a
type PagReqHeaders ranges a = Header "Range" (AppliedRanges ranges a)
type PagRespHeaders ranges a = Header "Total-Count" Int ': PageHeaders ranges a

type DisruptionFields = '["id", "occurredOn", "updatedOn", "resolvedOn", "duration"]
type DisruptionRanges = AppliedRanges DisruptionFields UIDisruption
type DisruptionPagReqHeaders = PagReqHeaders DisruptionFields UIDisruption
type DisruptionPagRespHeaders = PagRespHeaders DisruptionFields UIDisruption

type DataApi =
       "stats" :> StatsApi
  :<|> "disruptions" :> DisruptionApi
  :<|> "objects" :> ObjectsApi
  :<|> "facilities" :> FacilitiesApi

type StatsApi =
       Get '[JSON] UIOverallStats
  :<|> "averageDisruptionsPerDay" :> Get '[JSON] [UIDisruptionPerDay]

type DisruptionApi =
       DisruptionPagReqHeaders :> GetPartialContent '[JSON] (Headers DisruptionPagRespHeaders [UIDisruption])
  :<|> "markers" :> Get '[JSON] [UIMapMarker]
  :<|> "active" :> DisruptionPagReqHeaders :> GetPartialContent '[JSON] (Headers DisruptionPagRespHeaders [UIDisruption])

type ObjectsApi = QueryParam "search" Text :> Get '[JSON] [UIObjectSearchResult]

type FacilitiesApi =
       Capture "facilityId" FacilityId :> Get '[JSON] UIFacilityDetails
  :<|> DisruptionPagReqHeaders :> Capture "facilityId" FacilityId :> "disruptions" :>  GetPartialContent '[JSON] (Headers DisruptionPagRespHeaders [UIDisruption])

dataServer ::
     IORef UIOverallStats
  -> IORef DisruptionsPerDay
  -> IORef [UIDisruption]
  -> IORef [UIMapMarker]
  -> IORef [UIObjectSearchResult]
  -> IORef (Map FacilityId UIFacilityDetails)
  -> Server DataApi
dataServer stats disruptionsPerDay disruptions marker searchResults facilities =
  statsApi stats disruptionsPerDay :<|> disruptionsApi disruptions marker :<|> objectsApi searchResults :<|> facilitiesApi facilities disruptions

statsApi :: IORef UIOverallStats -> IORef DisruptionsPerDay -> Server StatsApi
statsApi statsRef disPerDayRef =
  readIORef statsRef :<|> disPerDayHandler
  where
    disPerDayHandler = reverse . fmap mkUIDisruptionPerDay . take 31 . drop 1 . reverse . M.toList . fmap average .  snd <$> readIORef disPerDayRef

disruptionsApi :: IORef [UIDisruption] -> IORef [UIMapMarker] -> Server DisruptionApi
disruptionsApi disRef markerRef =
  allDisruptionsHandler :<|> readIORef markerRef :<|> activeDisruptionsHandler
  where
    allDisruptionsHandler ranges = readIORef disRef >>= disruptionsHandler ranges
    activeDisruptionsHandler ranges = readIORef disRef >>= disruptionsHandler ranges . filter (isNothing . uidResolvedOn)

objectsApi :: IORef [UIObjectSearchResult] -> Server ObjectsApi
objectsApi ref maybeSearch =
  if length search < 3
  then return []
  else filter (nameMatches . uioName . uiosrObject) <$> liftIO (readIORef ref)

  where
    search = T.toLower $ fromMaybe "" maybeSearch
    nameMatches = T.isInfixOf search . T.toLower

facilitiesApi :: IORef (Map FacilityId UIFacilityDetails) -> IORef [UIDisruption] -> Server FacilitiesApi
facilitiesApi facilities disruptions =
  facilityDetailsHandler :<|> facilityDisruptionsHandler
  where
    facilityDetailsHandler facilityId =
      maybe (throwError err404) return . lookup facilityId =<< readIORef facilities
    facilityDisruptionsHandler ranges facilityId =
      readIORef disruptions >>= disruptionsHandler ranges . filter ((==) facilityId . uidFacilityId)

disruptionsHandler :: Maybe DisruptionRanges -> [UIDisruption] -> Handler (Headers DisruptionPagRespHeaders [UIDisruption])
disruptionsHandler r ds =
  addHeader (length ds) <$>
    fromMaybe (returnIdRange ds defaultRange)
        (fmap (returnIdRange ds) (r >>= extractRange)
    <|> fmap (returnOccurredOnRange ds) (r >>= extractRange)
    <|> fmap (returnUpdatedOnRange ds) (r >>= extractRange)
    <|> fmap (returnResolvedOnRange ds) (r >>= extractRange)
    <|> fmap (returnDurationRange ds) (r >>= extractRange))

  where
    returnIdRange xs (range :: Range "id" Int) =
      returnRange range (applyRange range xs)

    returnOccurredOnRange xs (range :: Range "occurredOn" DateTime) =
      returnRange range (applyRange range xs)

    returnUpdatedOnRange xs (range :: Range "updatedOn" (Maybe DateTime)) =
      returnRange range (applyRange range xs)

    returnResolvedOnRange xs (range :: Range "resolvedOn" (Maybe DateTime)) =
      returnRange range (applyRange range xs)

    returnDurationRange xs (range :: Range "duration" Integer) =
      returnRange range (applyRange range xs)

defaultRange :: Range "id" Int
defaultRange = getDefaultRange (Proxy @UIDisruption)

max_stars_repo_name: Akii/elescore

max_issues_repo_name: Akii/elescore

max_forks_repo_name: Akii/elescore