text:
-------------------------------------------------------------------------------
-- |
-- Module    :  Classifiers.RL.Control.MonadMWCRandom
-- Copyright :  (c) Sentenai 2017
-- License   :  BSD3
-- Maintainer:  sam@sentenai.com
-- Stability :  experimental
-- Portability: non-portable
--
-- typeclass to remove extraneous mwc-random functions
-------------------------------------------------------------------------------
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE UndecidableInstances #-}
module Control.MonadMWCRandom
  ( MonadMWCRandom(..)
  , MonadMWCRandomIO
  , MWCRand
  , MWCRandT(..)
  , runMWCRand
  , runMWCRandT
  -- * re-exports from System.Random.MWC
  , GenIO
  -- * wrappers for System.Random.MWC
  , uniform
  , uniformR
  -- * wrappers for Statistics.Distribution
  , Control.MonadMWCRandom.genContVar
  -- * extras
  , sampleFrom
  , Variate
  , _uniform
  ) where

import System.Random.MWC (GenIO, Variate)
import qualified System.Random.MWC as MWC
import qualified Statistics.Distribution as Stats
import Control.MonadEnv (MonadEnv(..), Obs, Initial)
import Control.Monad.Identity
import Control.Monad.IO.Class
import Control.Monad.Trans.Reader
import Control.Monad.Trans.Writer
import Control.Monad.Trans.State
import Control.Monad.Trans
import Control.Monad.Trans.RWS (RWST)
import Control.Exception.Safe
import Control.Monad.Primitive (PrimState, PrimMonad)

-- | a convenience helper to reference the underlying System.Random.MWC function
_uniform :: (PrimMonad m, Variate a) => MWC.Gen (PrimState m) -> m a
_uniform = MWC.uniform

-- | MonadMWCRandom for public use. FIXME: use with PrimState so that we can use ST
class Monad m => MonadMWCRandom m where
  getGen :: m GenIO

-- | A convenience type constraint with MonadMWCRandom and MonadIO.
type MonadMWCRandomIO m = (MonadIO m, MonadMWCRandom m)

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

instance MonadMWCRandom m => MonadMWCRandom (StateT s m) where
  getGen :: StateT s m GenIO
  getGen = lift getGen

instance MonadMWCRandom m => MonadMWCRandom (ReaderT s m) where
  getGen :: ReaderT s m GenIO
  getGen = lift getGen

instance (Monoid w, MonadMWCRandom m) => MonadMWCRandom (WriterT w m) where
  getGen :: WriterT w m GenIO
  getGen = lift getGen

instance (Monoid w, MonadMWCRandom m) => MonadMWCRandom (RWST r w s m) where
  getGen :: RWST r w s m GenIO
  getGen = lift getGen

-- | in the end, we can always use IO to get our generator, but we will create a
-- new generator on each use.
instance MonadMWCRandom IO where
  getGen :: IO GenIO
  getGen = MWC.createSystemRandom

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

-- | uniform referencing MonadMWCRandom's generator
uniform :: (MonadIO m, MonadMWCRandom m, Variate a) => m a
uniform = getGen >>= liftIO . MWC.uniform

-- | uniformR referencing MonadMWCRandom's generator
uniformR :: (MonadIO m, MonadMWCRandom m, Variate a) => (a, a) -> m a
uniformR r = getGen >>= liftIO . MWC.uniformR r

-- | genContVar referencing MonadMWCRandom's generator
genContVar :: (MonadIO m, MonadMWCRandom m, Stats.ContGen d) => d -> m Double
genContVar d = getGen >>= liftIO . Stats.genContVar d

-- ========================================================================= --
-- * Utility functions functions

-- | Sample a single index from a list of weights, converting the list into
-- a distribution
sampleFrom :: (MonadIO m, MonadMWCRandom m) => [Double] -> m (Int, [Double])
sampleFrom xs = fmap ((,dist) . choose) uniform
  where
    dist :: [Double]
    dist = fmap (/ total) xs

    total :: Double
    total = sum xs

    choose :: Double -> Int
    choose n =
      -- Return the head index (unsafeHead is safe since the last elem's snd must be 1.0)
      fst . head .

      -- Drop while the cumulative sum is < the given value
      dropWhile ((< n) . snd) .

      -- Pair each elem with its index
      zip [0..] .

      -- Transform list of probabilities to cumulative sum
      scanl1 (+) $ dist

-- ========================================================================= --
-- * A concrete type for MonadMWCRandom

-- | a wrapper to share a generator without using reader
newtype MWCRandT m a = MWCRandT { getMWCRandT :: ReaderT GenIO m a }
  deriving (Functor, Applicative, Monad, MonadTrans, MonadThrow, MonadIO)

-- | unwrap MonadMWCRandom
runMWCRandT :: MWCRandT m a -> GenIO -> m a
runMWCRandT = runReaderT . getMWCRandT

-- | simple type alias for transformer-less variant
type MWCRand = MWCRandT Identity

-- | run a transformerless MWC-random Monad
runMWCRand :: MWCRand a -> GenIO -> a
runMWCRand = runMWCRand

-- | instance declaration of MonadMWCRandom for MWCRandT
instance Monad m => MonadMWCRandom (MWCRandT m) where
  getGen :: MWCRandT m GenIO
  getGen = MWCRandT ask

-- | An instance which allows for an environment to hold a reference to a shared
-- MWC-random generator
instance MonadEnv m s a r => MonadEnv (MWCRandT m) s a r where
  reset :: MWCRandT m (Initial s)
  reset = lift reset

  step :: a -> MWCRandT m (Obs r s)
  step a = lift $ step a

max_stars_repo_name: juliendehos/reinforce

max_issues_repo_name: sentenai/reinforce

max_forks_repo_name: sentenai/reinforce


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Control.DeepSeq
-- Copyright   :  (c) The University of Glasgow 2001-2009
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  libraries@haskell.org
-- Stability   :  stable
-- Portability :  portable
--
-- Provides an overloaded function 'deepseq' for fully evaluating data
-- structures.
--

module Control.DeepSeq (
     deepseq,
     NFData(..),
  ) where

import Data.Int
import Data.Word
import Data.Ratio
import Data.Complex
import Data.Map
import Data.Set
import Data.IntMap
import Data.IntSet
import Data.Tree
import Data.Array

-- | Fully evaluates its argument.  The name 'deepseq' is used to
-- illustrate the relationship to 'seq': where 'seq' is shallow in
-- the sense that it only evaluates the top level of its argument,
-- 'deepseq' traverses the entire data structure evaluating it
-- completely.
--
-- 'deepseq' can be useful for forcing pending exceptions,
-- eradicating space leaks, or forcing lazy I/O to happen.  It is
-- also useful in conjunction with parallel Strategies (see the
-- @parallel@ package).
--
-- There is no guarantee about the ordering of evaluation.  The
-- implementation may evaluate the components of the structure in
-- any order or in parallel.  To impose an actual order on
-- evaluation, use 'pseq' from "Control.Parallel" in the
-- @parallel@ package.
--
deepseq :: NFData a => a -> b -> b
deepseq a b = rnf a `seq` b

-- A class of types that can be fully evaluated.
class NFData a where
    -- | rnf should reduce its argument to normal form (that is, fully
    -- evaluate all sub-components), and then return '()'.
    --
    -- The default implementation of 'rnf' is
    --
    -- > rnf a = a `seq` ()
    --
    -- which may be convenient when defining instances for data types with
    -- no unevaluated fields (e.g. enumerations).
    rnf :: a -> ()
    rnf a = a `seq` ()

instance NFData Int
instance NFData Word
instance NFData Integer
instance NFData Float
instance NFData Double

instance NFData Char
instance NFData Bool
instance NFData ()

instance NFData Int8
instance NFData Int16
instance NFData Int32
instance NFData Int64

instance NFData Word8
instance NFData Word16
instance NFData Word32
instance NFData Word64

--Rational and complex numbers.

instance (Integral a, NFData a) => NFData (Ratio a) where
  rnf x = rnf (numerator x, denominator x)

instance (RealFloat a, NFData a) => NFData (Complex a) where
  rnf (x:+y) = rnf x `seq`
               rnf y `seq`
               ()

instance NFData a => NFData (Maybe a) where
    rnf Nothing  = ()
    rnf (Just x) = rnf x

instance (NFData a, NFData b) => NFData (Either a b) where
    rnf (Left x)  = rnf x
    rnf (Right y) = rnf y

instance (NFData k, NFData a) => NFData (Data.Map.Map k a) where
    rnf = rnf . Data.Map.toList

instance NFData a => NFData (Data.Set.Set a) where
    rnf = rnf . Data.Set.toList

instance NFData a => NFData (Data.Tree.Tree a) where
    rnf (Data.Tree.Node r f) = rnf r `seq` rnf f

instance NFData a => NFData (Data.IntMap.IntMap a) where
    rnf = rnf . Data.IntMap.toList

instance NFData Data.IntSet.IntSet where
    rnf = rnf . Data.IntSet.toList

instance NFData a => NFData [a] where
    rnf [] = ()
    rnf (x:xs) = rnf x `seq` rnf xs

instance (Ix a, NFData a, NFData b) => NFData (Array a b) where
    rnf x = rnf (bounds x, Data.Array.elems x)

instance (NFData a, NFData b) => NFData (a,b) where
  rnf (x,y) = rnf x `seq` rnf y

instance (NFData a, NFData b, NFData c) => NFData (a,b,c) where
  rnf (x,y,z) = rnf x `seq` rnf y `seq` rnf z

instance (NFData a, NFData b, NFData c, NFData d) => NFData (a,b,c,d) where
  rnf (x1,x2,x3,x4) = rnf x1 `seq`
		        rnf x2 `seq`
		        rnf x3 `seq`
		        rnf x4

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5) =>
         NFData (a1, a2, a3, a4, a5) where
  rnf (x1, x2, x3, x4, x5) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6) =>
         NFData (a1, a2, a3, a4, a5, a6) where
  rnf (x1, x2, x3, x4, x5, x6) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6, NFData a7) =>
         NFData (a1, a2, a3, a4, a5, a6, a7) where
  rnf (x1, x2, x3, x4, x5, x6, x7) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6 `seq`
                  rnf x7

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6, NFData a7, NFData a8) =>
         NFData (a1, a2, a3, a4, a5, a6, a7, a8) where
  rnf (x1, x2, x3, x4, x5, x6, x7, x8) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6 `seq`
                  rnf x7 `seq`
                  rnf x8

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6, NFData a7, NFData a8, NFData a9) =>
         NFData (a1, a2, a3, a4, a5, a6, a7, a8, a9) where
  rnf (x1, x2, x3, x4, x5, x6, x7, x8, x9) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6 `seq`
                  rnf x7 `seq`
                  rnf x8 `seq`
                  rnf x9

ext: hs


text:
module Main where

import Plot
import Graphics.GD
import Data.Complex

imageDim :: Size
imageDim = (1200, 400)
-- 7200, 4800

--imageDim :: Size
--imageDim = (600, 600)

-- How many iterations before  we declare the point doesn't diverge
maxIter = 64

-- Coordinates in the complex plane to graph
-- First coordinate must be top left.
mbWindow :: Window
mbWindow = ( (-3,1), (3,-1) )
-- mbWindow = ( (-3,2), (3, -2))

--mbWindow :: Window
--mbWindow = ( (-1.5,1.5), (1.5,-1.5) ) -- Full view

-- C value
c :: Complex Double
c = (-0.423) :+ 0.745 -- Dust
-- c = (0.45) :+ (-0.1428)
-- c = (0.285 :+ 0.01)
-- c = (-0.45) :+ (-0.143) -- hazycloud
-- c = (-0.54) :+ (0.54) -- dragon
-- c = 0 :+ (0.8) -- veins
-- c = 0.355 :+ (0.355) -- swirlybois

main :: IO ()
main = drawPlot drawJulia imageDim mbWindow "juliaset.png"

drawJulia :: Coordinate -> Color
drawJulia (x,y) = colorIterations $ julia c (x :+ y) 0

-- Count number of iterations for the given point to diverge
-- Non-diverging points return 0
julia :: Complex Double -- Constant C
      -> Complex Double -- Iterating Z value
      -> Int -- Current iteration
      -> Int -- Iterations before diverging
julia c z iter
    | iter > maxIter = 0
    | otherwise = let z' = z^2 + c in
                  if magnitude z' > 2
                  then iter
                  else julia c z' (iter+1)

-- original colors
--colorIterations :: Int -> Color
--colorIterations x
--    | x > maxIter = rgb 255 255 255
--    | otherwise = let c = x*4
--                  in rgb x c x

-- the version we used for the purple julias
colorIterations :: Int -> Color
colorIterations x
    | x > maxIter = rgb 255 255 255
    | otherwise = let c = 50 + x
                      x'' = 50 + x*3 -- addition because to avoid fading to black
                      x' = 100 + x*2
                  in rgb x'' c x'

max_stars_repo_name: typeclasses/fractal-sets

max_issues_repo_name: typeclasses/fractal-sets

max_forks_repo_name: typeclasses/fractal-sets


text:
{-# LANGUAGE DataKinds, TypeFamilies #-}
{-# OPTIONS_GHC -Wno-missing-export-lists #-}
-- | Scalar-based implementation of fully connected neutral network
-- for classification of MNIST digits. Sports 2 hidden layers.
module HordeAd.Tool.MnistFcnnScalar where

import Prelude

import           Control.Exception (assert)
import           Data.Proxy (Proxy)
import qualified Data.Strict.Vector as Data.Vector
import qualified Data.Vector.Generic as V
import           GHC.Exts (inline)
import           Numeric.LinearAlgebra (Vector)

import HordeAd.Core.DualNumber
import HordeAd.Core.Engine
import HordeAd.Core.PairOfVectors (DualNumberVariables, var0)
import HordeAd.Tool.MnistData

-- | Compute the output of a neuron, without applying activation function,
-- from trainable inputs in @xs@ and parameters (the bias and weights)
-- at @variables@ starting at @offset@. Useful for neurons in the middle
-- of the network, receiving inputs from other neurons.
sumTrainableInputs
  :: forall d r m. DualMonad d r m
  => Data.Vector.Vector (DualNumber d r) -> Int -> DualNumberVariables d r
  -> m (DualNumber d r)
sumTrainableInputs xs offset variables = do
  let bias = var0 variables offset
      f :: DualNumber d r -> Int -> DualNumber d r -> DualNumber d r
      f !acc i u =
        let v = var0 variables (offset + 1 + i)
        in acc + u * v
  returnLet $ V.ifoldl' f bias xs

-- | Compute the output of a neuron, without applying activation function,
-- from constant data in @xs@ and parameters (the bias and weights)
-- at @variables@ starting at @offset@. Useful for neurons at the bottom
-- of the network, tasked with ingesting the data.
sumConstantData
  :: forall d r m. DualMonad d r m
  => Vector r -> Int -> DualNumberVariables d r -> m (DualNumber d r)
sumConstantData xs offset variables = do
  let bias = var0 variables offset
      f :: DualNumber d r -> Int -> r -> DualNumber d r
      f !acc i r =
        let v = var0 variables (offset + 1 + i)
        in acc + scale r v
  returnLet $ V.ifoldl' f bias xs

hiddenLayerMnist
  :: forall d r m. DualMonad d r m
  => (DualNumber d r -> m (DualNumber d r)) -> Vector r
  -> DualNumberVariables d r -> Int
  -> m (Data.Vector.Vector (DualNumber d r))
hiddenLayerMnist factivation input variables width = do
  let nWeightsAndBias = V.length input + 1
      f :: Int -> m (DualNumber d r)
      f i = do
        outSum <- sumConstantData input (i * nWeightsAndBias) variables
        factivation outSum
  V.generateM width f

middleLayerMnist
  :: forall d r m. DualMonad d r m
  => (DualNumber d r -> m (DualNumber d r))
  -> Data.Vector.Vector (DualNumber d r)
  -> Int -> DualNumberVariables d r -> Int
  -> m (Data.Vector.Vector (DualNumber d r))
middleLayerMnist factivation hiddenVec offset variables width = do
  let nWeightsAndBias = V.length hiddenVec + 1
      f :: Int -> m (DualNumber d r)
      f i = do
        outSum <- sumTrainableInputs hiddenVec
                                     (offset + i * nWeightsAndBias)
                                     variables
        factivation outSum
  V.generateM width f

outputLayerMnist
  :: forall d r m. DualMonad d r m
  => (Data.Vector.Vector (DualNumber d r)
      -> m (Data.Vector.Vector (DualNumber d r)))
  -> Data.Vector.Vector (DualNumber d r) -> Int
  -> DualNumberVariables d r -> Int
  -> m (Data.Vector.Vector (DualNumber d r))
outputLayerMnist factivation hiddenVec offset variables width = do
  let nWeightsAndBias = V.length hiddenVec + 1
      f :: Int -> m (DualNumber d r)
      f i = sumTrainableInputs hiddenVec
                               (offset + i * nWeightsAndBias)
                               variables
  vOfSums <- V.generateM width f
  factivation vOfSums

fcnnMnistLen0 :: Int -> Int -> Int
fcnnMnistLen0 widthHidden widthHidden2 =
  widthHidden * (sizeMnistGlyph + 1)
  + widthHidden2 * (widthHidden + 1)
  + sizeMnistLabel * (widthHidden2 + 1)

-- | Fully connected neural network for the MNIST digit classification task.
-- There are two hidden layers and both use the same activation function.
-- The output layer uses a different activation function.
-- The widths of the hidden layers are @widthHidden@ and @widthHidden2@
-- and from these, the @fcnnMnistLen2@ function computes the number
-- of scalar dual number parameters (variables) to be given to the program.
fcnnMnist0 :: forall d r m. DualMonad d r m
           => (DualNumber d r -> m (DualNumber d r))
           -> (Data.Vector.Vector (DualNumber d r)
               -> m (Data.Vector.Vector (DualNumber d r)))
           -> Int
           -> Int
           -> Vector r
           -> DualNumberVariables d r
           -> m (Data.Vector.Vector (DualNumber d r))
fcnnMnist0 factivationHidden factivationOutput widthHidden widthHidden2
           input variables = do
  let !_A = assert (sizeMnistGlyph == V.length input) ()
  layer1 <- inline hiddenLayerMnist factivationHidden input
                                    variables widthHidden
  let offsetMiddle = widthHidden * (sizeMnistGlyph + 1)
  layer2 <- inline middleLayerMnist factivationHidden layer1
                                    offsetMiddle variables widthHidden2
  let offsetOutput = offsetMiddle + widthHidden2 * (widthHidden + 1)
  inline outputLayerMnist factivationOutput layer2
                          offsetOutput variables sizeMnistLabel

-- | The neural network applied to concrete activation functions
-- and composed with the appropriate loss function.
fcnnMnistLoss0
  :: DualMonad d r m
  => Int -> Int -> MnistData r -> DualNumberVariables d r
  -> m (DualNumber d r)
fcnnMnistLoss0 widthHidden widthHidden2 (input, target) variables = do
  result <- inline fcnnMnist0 logisticAct softMaxAct
                              widthHidden widthHidden2 input variables
  lossCrossEntropy target result

-- | A function testing the neural network given testing set of inputs
-- and the trained parameters.
--
-- The proxy argument is needed only for the (spurious) SPECIALIZE pragma,
-- becuase I can't write @SPECIALIZE fcnnMnistTest0 \@Double@.
fcnnMnistTest0 :: forall r. IsScalar 'DModeGradient r
           => Proxy r -> Int -> Int -> [MnistData r] -> Domain0 r
           -> r
fcnnMnistTest0 _ widthHidden widthHidden2 inputs params0 =
  let matchesLabels :: MnistData r -> Bool
      matchesLabels (glyph, label) =
        let nn = inline (fcnnMnist0 @'DModeGradient) logisticAct softMaxAct
                                        widthHidden widthHidden2 glyph
            value = V.map (\(D r _) -> r)
                    $ primalValueGeneral nn (params0, V.empty, V.empty, V.empty)
        in V.maxIndex value == V.maxIndex label
  in fromIntegral (length (filter matchesLabels inputs))
     / fromIntegral (length inputs)
{-# SPECIALIZE fcnnMnistTest0 :: Proxy Double -> Int -> Int -> [MnistData Double] -> Domain0 Double -> Double #-}

max_stars_repo_name: Mikolaj/horde-ad

max_issues_repo_name: Mikolaj/horde-ad

max_forks_repo_name: Mikolaj/horde-ad


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RecordWildCards #-}
module NanoML.Learn where

import Control.Arrow (second)
import Control.Monad
import Data.Aeson
import qualified Data.ByteString.Lazy as LBS
import qualified Data.Csv as CSV
import Data.Function
import Data.List
import qualified Data.List as List
import Data.Maybe
import qualified Data.Set as Set
import qualified Data.Vector as V
import GHC.Generics
-- import Grenade
import Numeric.LinearAlgebra hiding (rank)
import Text.Read

import NanoML.Classify
import NanoML.Monad
import NanoML.Types

rankExprs :: Net -> [Feature] -> Prog -> [(Confidence, Constraint)]
rankExprs net fs p = second getThing <$!> rank net samples
  where
  samples = [ s {getThing = c}
            | s <- concatMap mkfsD tbad
            , c <- maybeToList (List.find (\c -> getThing s == fromJust (constraintSpan c)) cores)
            ]

  (tbad, cores, me) = case runEval stdOpts (typeProg p) of
    Left e -> ([], [], Just e)
    Right (tp, cs) -> (tp, Set.toList (mconcat cs), Nothing)

  mkfsD (TDFun _ _ pes) = mconcat (map (mkTypeOut . snd) pes)
  mkfsD (TDEvl _ e)     = mkTypeOut e
  mkfsD _               = mempty

  -- inSlice s = any (\c -> getSpan s == fromJust (constraintSpan c)) cores

  -- mkTypeOut :: TExpr -> [Sample ()]
  mkTypeOut te = actfold f [] te
    where
    f p e acc = (:acc) $ MkSample
                { getThing = infoSpan (texprInfo e)
                , getSample = vector $ concatMap (\(_,c) -> c p e) fs
                }

stdNet :: IO (Net, [Feature])
stdNet = loadNet "models/op+context+type-hidden-500.json" >>= \case
  Nothing -> fail "could not decode 'models/op+context+type-hidden-500.json'"
  Just net -> return (net, preds_tis ++ map only_ctx preds_tis_ctx ++ preds_tcon_ctx)

data Prediction = Change | NoChange
  deriving (Eq, Show)

-- predict :: Net -> Sample -> Prediction
-- predict = undefined

type Confidence = Double

rank :: Net -> [Sample a] -> [(Confidence, Sample a)]
rank net samples =
  sortBy (flip compare `on` fst)
    [ (r ! 1, s)
    | s <- samples
    , let r = runNet net (getSample s)
    ]

data Net = MkNet
  { hidden :: [Weights]
  , output :: Weights
  } deriving (Eq, Show, Generic)
instance FromJSON Net

data Weights = MkWeights
  { weights :: Matrix Double
  , biases :: Vector Double
  } deriving (Eq, Show, Generic)
instance FromJSON Weights where
  parseJSON = withObject "Weights" $ \v -> MkWeights
    <$> fmap mkmatrix (v .: "weights")
    <*> fmap vector (v .: "biases")
    where
    mkmatrix (ws :: [[Double]]) = matrix (length $ head ws) (concat ws)

loadNet :: FilePath -> IO (Maybe Net)
loadNet file = decode <$> LBS.readFile file

runNet :: Net -> Vector Double -> Vector Double
runNet MkNet{..} features =
  softmax (runLayer (foldl' runHidden features hidden) output)
  where
  runHidden fs ws = relu $ runLayer fs ws
  runLayer fs MkWeights{..} =
    fs <# weights + biases

relu :: Vector Double -> Vector Double
relu v = cmap (\x -> max 0 x) v

softmax :: Vector Double -> Vector Double
softmax v = cmap (\x -> exp x / summed) v
  where
  summed = sumElements (cmap exp v)

data Sample a = MkSample { getThing :: a, getSample :: Vector Double }
  deriving (Eq, Show, Generic)
instance CSV.FromRecord (Sample SrcSpan) where
  parseRecord r = MkSample
    <$> (maybe mzero return . readMaybe =<< CSV.parseField (r V.! 0))
    <*> fmap vector (CSV.parseRecord (V.drop 4 r))

loadSamples :: FilePath -> IO (Either String (V.Vector (Sample SrcSpan)))
loadSamples file = CSV.decode CSV.HasHeader <$> LBS.readFile file

-- net2network :: Net -> Network ts ds
-- net2network MkNet{..} = foldr addLayer NNil (output : reverse hidden)

-- addLayer :: Weights -> Network ts ds
--          -> Network (FullyConnected dx dy : Relu : ts) (d : d : ds)
-- addLayer ws net = mkLayer ws :~> Relu :~> net

max_stars_repo_name: gsakkas/nanomaly

max_issues_repo_name: gsakkas/nanomaly

max_forks_repo_name: gsakkas/nanomaly


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts  #-}
{-# LANGUAGE ViewPatterns      #-}
module Tests.Transform
    (
      tests
    ) where

import Data.Bits ((.&.), shiftL)
import Data.Complex (Complex((:+)))
import Data.Functor ((<$>))
import Numeric.Sum (kbn, sumVector)
import Statistics.Function (within)
import Statistics.Transform (CD, dct, fft, idct, ifft)
import Test.Framework (Test, testGroup)
import Test.Framework.Providers.QuickCheck2 (testProperty)
import Test.QuickCheck ( Positive(..), Arbitrary(..), Blind(..), (==>), Gen
                       , choose, vectorOf, counterexample, forAll)
import Test.QuickCheck.Property (Property(..))
import Tests.Helpers (testAssertion)
import Text.Printf (printf)
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U

tests :: Test
tests = testGroup "fft" [
          testProperty "t_impulse"        t_impulse
        , testProperty "t_impulse_offset" t_impulse_offset
        , testProperty "ifft . fft = id"  (t_fftInverse $ ifft . fft)
        , testProperty "fft . ifft = id"  (t_fftInverse $ fft . ifft)
        , testProperty "idct . dct = id [up to scale]"
            (t_fftInverse (\v -> U.map (/ (2 * fromIntegral (U.length v))) $ idct $ dct v))
        , testProperty "dct . idct = id [up to scale]"
            (t_fftInverse (\v -> U.map (/ (2 * fromIntegral (U.length v))) $ idct $ dct v))
          -- Exact small size DCT
          -- 1
        , testDCT [1] $ [2]
          -- 2
        , testDCT [1,0] $ map (*2) [1, cos (pi/4)   ]
        , testDCT [0,1] $ map (*2) [1, cos (3*pi/4) ]
          -- 4
        , testDCT [1,0,0,0] $ map (*2) [1, cos(  pi/8), cos( 2*pi/8), cos( 3*pi/8)]
        , testDCT [0,1,0,0] $ map (*2) [1, cos(3*pi/8), cos( 6*pi/8), cos( 9*pi/8)]
        , testDCT [0,0,1,0] $ map (*2) [1, cos(5*pi/8), cos(10*pi/8), cos(15*pi/8)]
        , testDCT [0,0,0,1] $ map (*2) [1, cos(7*pi/8), cos(14*pi/8), cos(21*pi/8)]
          -- Exact small size IDCT
          -- 1
        , testIDCT [1] [1]
          -- 2
        , testIDCT [1,0]            [1,         1          ]
        , testIDCT [0,1] $ map (*2) [cos(pi/4), cos(3*pi/4)]
          -- 4
        , testIDCT [1,0,0,0]            [1,            1,            1,            1            ]
        , testIDCT [0,1,0,0] $ map (*2) [cos(   pi/8), cos( 3*pi/8), cos( 5*pi/8), cos( 7*pi/8) ]
        , testIDCT [0,0,1,0] $ map (*2) [cos( 2*pi/8), cos( 6*pi/8), cos(10*pi/8), cos(14*pi/8) ]
        , testIDCT [0,0,0,1] $ map (*2) [cos( 3*pi/8), cos( 9*pi/8), cos(15*pi/8), cos(21*pi/8) ]
        ]

-- A single real-valued impulse at the beginning of an otherwise zero
-- vector should be replicated in every real component of the result,
-- and all the imaginary components should be zero.
t_impulse :: Double -> Positive Int -> Bool
t_impulse k (Positive m) = U.all (c_near i) (fft v)
  where v = i `G.cons` G.replicate (n-1) 0
        i = k :+ 0
        n = 1 `shiftL` (m .&. 6)

-- If a real-valued impulse is offset from the beginning of an
-- otherwise zero vector, the sum-of-squares of each component of the
-- result should equal the square of the impulse.
t_impulse_offset :: Double -> Positive Int -> Positive Int -> Property
t_impulse_offset k (Positive x) (Positive m)
  -- For numbers smaller than 1e-162 their square underflows and test
  -- fails spuriously
  = abs k >= 1e-100 ==> U.all ok (fft v)
  where v = G.concat [G.replicate xn 0, G.singleton i, G.replicate (n-xn-1) 0]
        ok (re :+ im) = within ulps (re*re + im*im) (k*k)
        i  = k :+ 0
        xn = x `rem` n
        n  = 1 `shiftL` (m .&. 6)

-- Test that (ifft . fft ≈ id)
--
-- Approximate equality here is tricky. Smaller values of vector tend
-- to have large relative error. Thus we should test that vectors as
-- whole are approximate equal.
t_fftInverse :: (HasNorm (U.Vector a), U.Unbox a, Num a, Show a, Arbitrary a)
             => (U.Vector a -> U.Vector a) -> Property
t_fftInverse roundtrip =
  forAll (Blind <$> genFftVector) $ \(Blind x) ->
    let n  = G.length x
        x' = roundtrip x
        d  = G.zipWith (-) x x'
        nd = vectorNorm d
        nx = vectorNorm x
    in counterexample "Original vector"
     $ counterexample (show x )
     $ counterexample "Transformed one"
     $ counterexample (show x')
     $ counterexample (printf "Length = %i" n)
     $ counterexample (printf "|x - x'| / |x| = %.6g" (nd / nx))
     $ nd <= 3e-14 * nx

-- Test discrete cosine transform
testDCT :: [Double] -> [Double] -> Test
testDCT (U.fromList -> vec) (U.fromList -> res)
  = testAssertion ("DCT test for " ++ show vec)
  $ vecEqual 3e-14 (dct vec) res

-- Test inverse discrete cosine transform
testIDCT :: [Double] -> [Double] -> Test
testIDCT (U.fromList -> vec) (U.fromList -> res)
  = testAssertion ("IDCT test for " ++ show vec)
  $ vecEqual 3e-14 (idct vec) res

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

-- With an error tolerance of 8 ULPs, a million QuickCheck tests are
-- likely to all succeed. With a tolerance of 7, we fail around the
-- half million mark.
ulps :: Int
ulps = 8

c_near :: CD -> CD -> Bool
c_near (a :+ b) (c :+ d) = within ulps a c && within ulps b d

-- Arbitrary vector for FFT od DCT
genFftVector :: (U.Unbox a, Arbitrary a) => Gen (U.Vector a)
genFftVector = do
  n <- (2^)  <$> choose (1,9::Int)    -- Size of vector
  G.fromList <$> vectorOf n arbitrary -- Vector to transform

-- Ad-hoc type class for calculation of vector norm
class HasNorm a where
  vectorNorm :: a -> Double

instance HasNorm (U.Vector Double) where
  vectorNorm = sqrt . sumVector kbn . U.map (\x -> x*x)

instance HasNorm (U.Vector CD) where
  vectorNorm = sqrt . sumVector kbn . U.map (\(x :+ y) -> x*x + y*y)

-- Approximate equality for vectors
vecEqual :: Double -> U.Vector Double -> U.Vector Double -> Bool
vecEqual ε v u
  = vectorNorm (U.zipWith (-) v u) < ε * vectorNorm v

max_stars_repo_name: intricate/statistics

max_issues_repo_name: intricate/statistics

max_forks_repo_name: intricate/statistics


text:
module Polymorphic where

import Data.Complex
import Data.Ratio
import Prelude hiding (product)

stringLength :: String -> Int
stringLength x = case x of
                [] -> 0
                (x:xs) -> 1 + stringLength xs

integerListLength :: [Integer] -> Int
integerListLength x = case x of
                [] -> 0
                (x:xs) -> 1 + integerListLength xs

-- | A polymorphic length function
polymorphicLength :: [a] -> Int
polymorphicLength list = case list of
  []     -> 0
  _ : xs -> 1 + polymorphicLength xs

-- | A polymorphic reverse function
reverseOf :: [a] -> [a]
reverseOf xs = case xs of
            [] -> []
            c: cs -> reverseOf cs ++ [c]

-- | A polymorphic isPalindrome function
isPalindrome :: (Eq a) => [a]  -> Bool
isPalindrome xs
                     | null xs = True
                     | xs == reverseOf xs = True
                     | otherwise = False

-- | A polymorphic list equality function
listEqual :: (Eq a) => [a] -> [a] -> Bool
listEqual list1 list2 = case (list1, list2) of
  ([], [])       -> True
  ([], _ )       -> False
  (_ , [])       -> False
  (x: xs, y: ys) -> (x == y) && (xs `listEqual` ys)

isMonotonicallyIncreasing :: (Ord a) => [a]  -> Bool
isMonotonicallyIncreasing list = case list of
                             []         -> True
                             [x]        -> True
                             i1: i2: is | i1 <= i2  -> isMonotonicallyIncreasing (i2: is)
                                        | otherwise -> False

normaliseVector :: (Floating a) => [a] -> [a]
normaliseVector vector = divideByScalar vector (norm vector)

divideByScalar :: (Fractional a) => [a] -> a -> [a]
divideByScalar vector' scalar = case vector' of
                        []    -> []
                        f: fs -> (f / scalar): divideByScalar fs scalar

norm :: (Floating  a) => [a] -> a
norm vector' = sqrt (sumSqr vector')
     where
                 sumSqr :: (Num a) => [a] -> a
                 sumSqr vector'' = case vector'' of
                    []    -> 0
                    f: fs -> f*f + sumSqr fs

reallyPolymorphicLength :: (Integral b) => [a] -> b
reallyPolymorphicLength list = case list of
  []    -> 0
  _: xs -> 1 + reallyPolymorphicLength xs

printLargestAsString :: (Ord a, Show a) => [a] -> String
printLargestAsString list =
  "The largest element in the list is " ++ show (maximum list)

rationalZero :: Ratio Integer
floatZero    :: Float
doubleZero   :: Double
complexZero  :: Complex Double

rationalZero = 0
floatZero    = 0.0
doubleZero   = 0.0
complexZero  = mkPolar 0.0 0.0

product :: undefined -- TODO
product = undefined -- TODO


text:
module Scheme.Parser (readLisp, readLisps) where
import Control.Monad
import Control.Monad.Error (throwError)

import Data.Char
import Data.Maybe
import Data.Complex
import Data.Ratio
import Data.Vector as V (fromList)

import Text.Parsec hiding (spaces)

import Scheme.Types

readLisp :: FilePath -> String -> ThrowsErrorIO LispVal
readLisp fname inp = case parse (spaces >> lispVal) fname inp of
    Left err  -> throwError $ ParserError err
    Right val -> return val

readLisps :: FilePath -> String -> ThrowsErrorIO [LispVal]
readLisps fname inp = case parse lispParser fname inp of
    Left err  -> throwError $ ParserError err
    Right val -> return val

type Parser a = Parsec String () a

lispParser = spaces >> sepEndBy lispVal spaces

lispVal =
    try lispSequence <|> try num    <|> try lispChar <|> try bool <|>
    try atom         <|> stringLisp <|> quoted

comment = char ';' >> manyTill anyChar (try (char '\n' <|> char '\r'))

ignored = void space <|> void comment

spaces = many ignored

spaces1 :: Parser ()
spaces1 = skipMany1 ignored

symbol = oneOf "!$%&|*+-/:<=?>@^_~#"

bool = char '#' >> (string "t" <|> string "f")
                >>= (\s -> return $ case s of "t" -> Bool True
                                              "f" -> Bool False)

atom = do first <- letter <|> symbol
          rest  <- many (letter <|> symbol <|> digit)
          return . Atom $ first:rest

lispChar :: Parser LispVal
lispChar = do
    string "#\\"
    special <- optionMaybe $ choice (map (string . fst) specialCharLits)
    maybe (liftM Character anyChar) (\lit -> return (Character . fromJust $ lookup lit specialCharLits)) special
  where specialCharLits = [("newline", '\n'), ("space", ' '), ("tab", '\t'), ("return", '\r')]

num = do
    let radices = [('x', 16), ('d', 10), ('o', 8), ('t', 3), ('b', 2)]
    maybeBase <- optionMaybe (try (char '#' >> oneOf (map fst radices)))
    maybeSign <- optionMaybe (char '+' <|> char '-')
    let sign :: Num a => a -> a
        sign = fromMaybe id $ maybeSign >>= \c -> if c == '-' then Just negate else Nothing
        base = case maybeBase of
                    Nothing -> 10
                    Just s -> fromMaybe 10 (lookup s radices)
    real <- try (liftM (Double   . sign) (parseFloat base))
        <|> try (liftM (Rational . sign) (parseRational base))
        <|>      liftM (Int      . sign) (parseInteger base)
    maybeImaginary <- optionMaybe (try $ num >>= \n -> char 'i' >> return n)
    return $ maybe (Num real) (\(Num imag) ->
        Num $ Complex $ getAsDouble real :+ getAsDouble imag) maybeImaginary

digitCharacters = ['0'..'9'] ++ ['A'..'Z']

toDigit :: (Num a, Enum a) => Char -> a
toDigit c = case lookup (toUpper c) $ zip digitCharacters [0..] of
    Nothing -> error $ "character not between 0 and 9 or A and Z: '" ++ [c] ++ "'"
    Just d -> d

parseInteger :: Int -> Parser Integer
parseInteger base = liftM round $ parseDigits False base

parseFloat :: Int -> Parser Double
parseFloat base = do big <- parseDigits False base
                     char '.'
                     small <- parseDigits True base
                     return $ big + small

parseRational :: Int -> Parser (Ratio Integer)
parseRational base = do n <- parseInteger base
                        char '/'
                        d <- parseInteger base
                        return $ n % d

parseDigits :: (Floating a, Enum a) => Bool -> Int -> Parser a
parseDigits isDecimal base =
    let allowedDigits = take base digitCharacters
        convertStr numAsStr = snd $ foldr
            (\c (plc, val) -> (plc + 1, val + toDigit c * fromIntegral base ** plc))
            (if isDecimal then (-1) * fromIntegral (length numAsStr) else 0.0, 0.0) numAsStr
    in liftM convertStr $ many1 (satisfy (flip elem allowedDigits . toUpper))

escapables = [('t', '\t'), ('n', '\n'), ('r', '\r'), ('\\', '\\')]

stringLisp :: Parser LispVal
stringLisp = do char '"'
                s <- many (noneOf "\"\\" <|> (try $ do
                        char '\\'
                        esc <- anyChar
                        return $ fromMaybe esc $ lookup esc escapables))
                char '"'
                return $ String s

quoted = do char '\''
            v <- lispVal
            return $ List [Atom "quote", v]

lispSequence = do
    isVector <-  (string "("  >> return False)
             <|> (string "#(" >> return True)
    spaces
    vs <- sepEndBy lispVal spaces1
    maybeDot <- optionMaybe (spaces >> char '.' >> spaces1 >> lispVal)
    spaces
    char ')'
    return $ maybe
        (if isVector then Vector (toVector vs) else List vs)
        (DottedList vs) maybeDot

toVector = V.fromList

max_stars_repo_name: oroguh/SCMinHS

max_issues_repo_name: oroguh/SCMinHS

max_forks_repo_name: oroguh/SCMinHS


text:
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.FDistribution
-- Copyright : (c) 2011 Aleksey Khudyakov
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Fisher F distribution
module Statistics.Distribution.FDistribution (
    FDistribution
  , fDistribution
  , fDistributionNDF1
  , fDistributionNDF2
  ) where

import Data.Aeson (FromJSON, ToJSON)
import Data.Binary (Binary)
import Data.Data (Data, Typeable)
import Numeric.MathFunctions.Constants (m_neg_inf)
import GHC.Generics (Generic)
import qualified Statistics.Distribution as D
import Statistics.Function (square)
import Numeric.SpecFunctions (
  logBeta, incompleteBeta, invIncompleteBeta, digamma)
import Data.Binary (put, get)
import Control.Applicative ((<$>), (<*>))

-- | F distribution
data FDistribution = F { fDistributionNDF1 :: {-# UNPACK #-} !Double
                       , fDistributionNDF2 :: {-# UNPACK #-} !Double
                       , _pdfFactor        :: {-# UNPACK #-} !Double
                       }
                   deriving (Eq, Show, Read, Typeable, Data, Generic)

instance FromJSON FDistribution
instance ToJSON FDistribution

instance Binary FDistribution where
    get = F <$> get <*> get <*> get
    put (F x y z) = put x >> put y >> put z

fDistribution :: Int -> Int -> FDistribution
fDistribution n m
  | n > 0 && m > 0 =
    let n' = fromIntegral n
        m' = fromIntegral m
        f' = 0.5 * (log m' * m' + log n' * n') - logBeta (0.5*n') (0.5*m')
    in F n' m' f'
  | otherwise =
    error "Statistics.Distribution.FDistribution.fDistribution: non-positive number of degrees of freedom"

instance D.Distribution FDistribution where
  cumulative = cumulative

instance D.ContDistr FDistribution where
  density d x
    | x <= 0    = 0
    | otherwise = exp $ logDensity d x
  logDensity d x
    | x <= 0    = m_neg_inf
    | otherwise = logDensity d x
  quantile = quantile

cumulative :: FDistribution -> Double -> Double
cumulative (F n m _) x
  | x <= 0       = 0
  | isInfinite x = 1            -- Only matches +∞
  | otherwise    = let y = n*x in incompleteBeta (0.5 * n) (0.5 * m) (y / (m + y))

logDensity :: FDistribution -> Double -> Double
logDensity (F n m fac) x
  = fac + log x * (0.5 * n - 1) - log(m + n*x) * 0.5 * (n + m)

quantile :: FDistribution -> Double -> Double
quantile (F n m _) p
  | p >= 0 && p <= 1 =
    let x = invIncompleteBeta (0.5 * n) (0.5 * m) p
    in m * x / (n * (1 - x))
  | otherwise =
    error $ "Statistics.Distribution.Uniform.quantile: p must be in [0,1] range. Got: "++show p

instance D.MaybeMean FDistribution where
  maybeMean (F _ m _) | m > 2     = Just $ m / (m - 2)
                      | otherwise = Nothing

instance D.MaybeVariance FDistribution where
  maybeStdDev (F n m _)
    | m > 4     = Just $ 2 * square m * (m + n - 2) / (n * square (m - 2) * (m - 4))
    | otherwise = Nothing

instance D.Entropy FDistribution where
  entropy (F n m _) =
    let nHalf = 0.5 * n
        mHalf = 0.5 * m in
    log (n/m)
    + logBeta nHalf mHalf
    + (1 - nHalf) * digamma nHalf
    - (1 + mHalf) * digamma mHalf
    + (nHalf + mHalf) * digamma (nHalf + mHalf)

instance D.MaybeEntropy FDistribution where
  maybeEntropy = Just . D.entropy

instance D.ContGen FDistribution where
  genContVar = D.genContinous

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
module Kantour.DropCalc where

import System.Environment
import Control.Monad
import Control.Monad.Primitive
import System.Random.MWC
import Text.Printf
import qualified Data.Vector as V
import Statistics.Sample
import Data.Vector.Algorithms.Intro
import Data.Function
import Numeric.SpecFunctions
import Kantour.Subcommand

data SubCmdDropCalc

instance Subcommand SubCmdDropCalc where
    name _ = "DropCalc"
    main _ = defaultMain

experiment :: PrimMonad m => Gen (PrimState m) -> Double -> m Int
experiment g r = experiment' 0
  where
    experiment' cnt = do
        let newCnt = cnt+1
        v <- uniform g
        if v <= r
          then pure newCnt
          else experiment' newCnt

computeCDF :: Int -> Int -> Double -> Double
computeCDF k r p = i (1-p) (fromIntegral r) (fromIntegral k)
  where
    i x a b = incompleteBeta a b x

analyzeResult :: Int -> Double -> [Int] -> IO ()
analyzeResult _ _ [] = putStrLn "Result is empty."
analyzeResult tot p xs = do
    let _xs = V.fromList xs :: V.Vector Int
    ys' <- V.unsafeThaw _xs
    sort ys'
    -- now ys is the sorted version of xs
    ys <- V.unsafeFreeze ys'
    let xMin = V.head ys
        xMax = V.last ys
        ysDb = fromIntegral <$> ys :: V.Vector Double
    printf "min: %d, max: %d\n" xMin xMax
    printf "mean: %6.4f, std dev.: %6.4f\n" (mean ysDb) (stdDev ysDb)
    (fix $ \indScan curInd ->
      if curInd >= V.length ys
        then pure ()
        else do
          let curE = ys V.! curInd
          -- find next index whose elem is not curE
          let nextInd = fix (\self curInd' ->
                             if curInd' >= V.length ys
                               then Nothing
                               else if ys V.! curInd' /= curE
                                      then Just curInd'
                                      else self (curInd'+1)) (curInd+1)
              fI = fromIntegral :: Int -> Double
              ppr n cr = when (n `mod` 10 == 0) $
                  printf "<= %d\tactual: %6.4f%%\tcomputed: %6.4f%%\n"
                      n
                      (cr * 100)
                      (computeCDF n 1 (1-p) * 100)
          case nextInd of
              Nothing ->
                  ppr curE (fI (V.length ys -1) / fI tot )
              Just ind -> do
                  ppr curE (fI (ind-1) / fI tot)
                  indScan ind
          ) 0

defaultMain :: IO ()
defaultMain = do
    as <- getArgs
    case as of
        [rawDropRate,rawCnt]
            | Just (dropRate,cnt) <-
                (,) <$> parseRate rawDropRate
                    <*> parseCount rawCnt
            -> do
                printf "Performing %d experiments with drop rate %6.4f%%:\n"
                  cnt (dropRate * 100)
                g <- createSystemRandom
                results <- replicateM cnt (experiment g dropRate)
                analyzeResult cnt dropRate results
        _ -> putStrLn "Usage: <rate> <# of experiments>"
             >> putStrLn "example of rate: '0.20' or '20%'"
  where
    -- either a number between 0 and 1
    -- or {x}% where x can be between 0 and 100
    parseRate :: String -> Maybe Double
    parseRate [] = Nothing
    parseRate xs
        | last xs == '%'
        , [(r,"")] <- reads (init xs) = Just (r/100) >>= verify
        | [(r,"")] <- reads xs = Just r >>= verify
        | otherwise = Nothing
      where
        verify x = x <$ guard (x >= 0 && x <= 1)

    -- must be positive number
    parseCount :: String -> Maybe Int
    parseCount xs
        | [(v,"")] <- reads xs
        , v > 0 = Just v
        | otherwise = Nothing

max_stars_repo_name: Javran/tuppence

max_issues_repo_name: Javran/tuppence

max_forks_repo_name: Javran/tuppence


text:
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE Rank2Types            #-}
{-# LANGUAGE ViewPatterns          #-}
{-# OPTIONS_GHC -fno-warn-missing-methods #-}

-----------------------------------------------------------------------------
-- |
-- Module      :  Diagrams.TwoD.Apollonian
-- Copyright   :  (c) 2011, 2016 Brent Yorgey
-- License     :  BSD-style (see LICENSE)
-- Maintainer  :  byorgey@cis.upenn.edu
--
-- Generation of Apollonian gaskets.  Any three mutually tangent
-- circles uniquely determine exactly two others which are mutually
-- tangent to all three.  This process can be repeated, generating a
-- fractal circle packing.
--
-- See J. Lagarias, C. Mallows, and A. Wilks, \"Beyond the Descartes
-- circle theorem\", /Amer. Math. Monthly/ 109 (2002), 338--361.
-- <http://arxiv.org/abs/math/0101066>.
--
-- A few examples:
--
-- > import Diagrams.TwoD.Apollonian
-- > apollonian1 = apollonianGasket 0.01 2 2 2
--
-- <<diagrams/src_Diagrams_TwoD_Apollonian_apollonian1.svg#diagram=apollonian1&width=400>>
--
-- > import Diagrams.TwoD.Apollonian
-- > apollonian2 = apollonianGasket 0.01 2 3 3
--
-- <<diagrams/src_Diagrams_TwoD_Apollonian_apollonian2.svg#diagram=apollonian2&width=400>>
--
-- > import Diagrams.TwoD.Apollonian
-- > apollonian3 = apollonianGasket 0.01 2 4 7
--
-- <<diagrams/src_Diagrams_TwoD_Apollonian_apollonian3.svg#diagram=apollonian3&width=400>>
--
-----------------------------------------------------------------------------

module Diagrams.TwoD.Apollonian
       ( -- * Circles

         Circle(..), mkCircle, center, radius

         -- * Descartes' Theorem
       , descartes, other, initialConfig

         -- * Apollonian gasket generation

       , apollonian

         -- ** Kissing sets

       , KissingSet(..), kissingSets, flipSelected, selectOthers

         -- ** Apollonian trees

       , apollonianTrees, apollonianTree

         -- * Diagram generation

       , drawCircle
       , drawGasket
       , apollonianGasket

       ) where

import           Data.Complex
import qualified Data.Foldable    as F
import           Data.Maybe       (catMaybes)
import           Data.Tree

import           Diagrams.Prelude hiding (center, radius)

import           Control.Arrow    (second, (&&&))

------------------------------------------------------------
--  Circles
------------------------------------------------------------

-- | Representation for circles that lets us quickly compute an
--   Apollonian gasket.
data Circle n = Circle
  { bend :: n
    -- ^ The bend is the reciprocal of signed
    --   radius: a negative radius means the
    --   outside and inside of the circle are
    --   switched.  The bends of any four mutually
    --   tangent circles satisfy Descartes'
    --   Theorem.
  , cb   :: Complex n
    -- ^ /Product/ of bend and center represented
    --   as a complex number.  Amazingly, these
    --   products also satisfy the equation of
    --   Descartes' Theorem.
  }
  deriving (Eq, Show)

-- | Create a @Circle@ given a signed radius and a location for its center.
mkCircle :: Fractional n =>
            n -- ^ signed radius
         -> P2 n     -- ^ center
         -> Circle n
mkCircle r (unp2 -> (x,y)) = Circle (1/r) (b*x :+ b*y)
  where b = 1/r

-- | Get the center of a circle.
center :: Fractional n => Circle n -> P2 n
center (Circle b (cbx :+ cby)) = p2 (cbx / b, cby / b)

-- | Get the (unsigned) radius of a circle.
radius :: Fractional n => Circle n -> n
radius = abs . recip . bend

liftF :: RealFloat n => (forall a. Floating a => a -> a) -> Circle n -> Circle n
liftF f (Circle b c) = Circle (f b) (f c)

liftF2 :: RealFloat n => (forall a. Floating a => a -> a -> a) ->
          Circle n -> Circle n -> Circle n
liftF2 f (Circle b1 cb1) (Circle b2 cb2) = Circle (f b1 b2) (f cb1 cb2)

instance RealFloat n => Num (Circle n) where
  (+)           = liftF2 (+)
  (-)           = liftF2 (-)
  (*)           = liftF2 (*)
  negate        = liftF negate
  abs           = liftF abs
  fromInteger n = Circle (fromInteger n) (fromInteger n)

instance RealFloat n => Fractional (Circle n) where
  (/)   = liftF2 (/)
  recip = liftF recip

-- | The @Num@, @Fractional@, and @Floating@ instances for @Circle@
--   (all simply lifted elementwise over @Circle@'s fields) let us use
--   Descartes' Theorem directly on circles.
instance RealFloat n => Floating (Circle n) where
  sqrt = liftF sqrt

------------------------------------------------------------
--  Descartes' Theorem
------------------------------------------------------------

-- XXX generalize these for higher dimensions?

-- | Descartes' Theorem states that if @b1@, @b2@, @b3@ and @b4@ are
--   the bends of four mutually tangent circles, then
--
--   @
--     b1^2 + b2^2 + b3^2 + b4^2 = 1/2 * (b1 + b2 + b3 + b4)^2.
--   @
--
--   Surprisingly, if we replace each of the @bi@ with the /product/
--   of @bi@ and the center of the corresponding circle (represented
--   as a complex number), the equation continues to hold! (See the
--   paper referenced at the top of the module.)
--
--   @descartes [b1,b2,b3]@ solves for @b4@, returning both solutions.
--   Notably, @descartes@ works for any instance of @Floating@, which
--   includes both @Double@ (for bends), @Complex Double@ (for
--   bend/center product), and @Circle@ (for both at once).
descartes :: Floating n => [n] -> [n]
descartes [b1,b2,b3] = [r + s, -r + s]
  where r = 2 * sqrt (b1*b2 + b1*b3 + b2*b3)
        s = b1+b2+b3
descartes _ = error "descartes must be called on a list of length 3"

-- | If we have /four/ mutually tangent circles we can choose one of
--   them to replace; the remaining three determine exactly one other
--   circle which is mutually tangent.  However, in this situation
--   there is no need to apply 'descartes' again, since the two
--   solutions @b4@ and @b4'@ satisfy
--
--   @
--     b4 + b4' = 2 * (b1 + b2 + b3)
--   @
--
--   Hence, to replace @b4@ with its dual, we need only sum the other
--   three, multiply by two, and subtract @b4@.  Again, this works for
--   bends as well as bend/center products.
other :: Num n => [n] -> n -> n
other xs x = 2 * sum xs - x

-- | Generate an initial configuration of four mutually tangent
--   circles, given just the signed bends of three of them.
initialConfig :: RealFloat n => n -> n -> n -> [Circle n]
initialConfig b1 b2 b3 = cs ++ [c4]
  where cs     = [Circle b1 0, Circle b2 ((b2/b1 + 1) :+ 0), Circle b3 cb3]
        a      = 1/b1 + 1/b2
        b      = 1/b1 + 1/b3
        c      = 1/b2 + 1/b3
        x      = (b*b + a*a - c*c)/(2*a)
        y      = sqrt (b*b - x*x)
        cb3    = b3*x :+ b3*y
        [c4,_] = descartes cs

------------------------------------------------------------
--  Gasket generation
------------------------------------------------------------

select :: [a] -> [(a, [a])]
select [] = []
select (x:xs) = (x,xs) : (map . second) (x:) (select xs)

-- | The basic idea of a kissing set is supposed to represent a set of
--   four mutually tangent circles with one selected, though in fact
--   it is more general than that: it represents any set of objects
--   with one distinguished object selected.
data KissingSet n = KS { selected :: n, others :: [n] }
  deriving (Show)

-- | Generate all possible kissing sets from a set of objects by
--   selecting each object in turn.
kissingSets :: [n] -> [KissingSet n]
kissingSets = map (uncurry KS) . select

-- | \"Flip\" the selected circle to the 'other' circle mutually tangent
--   to the other three.  The new circle remains selected.
flipSelected :: Num n => KissingSet n -> KissingSet n
flipSelected (KS c cs) = KS (other cs c) cs

-- | Make the selected circle unselected, and select each of the
--   others, generating a new kissing set for each.
selectOthers :: KissingSet n -> [KissingSet n]
selectOthers (KS c cs) = [ KS c' (c:cs') | (c',cs') <- select cs ]

-- | Given a threshold radius and a list of /four/ mutually tangent
--   circles, generate the Apollonian gasket containing those circles.
--   Stop the recursion when encountering a circle with an (unsigned)
--   radius smaller than the threshold.
apollonian :: RealFloat n => n -> [Circle n] -> [Circle n]
apollonian thresh cs
  = (cs++)
  . concat
  . map (maybe [] flatten . prune p . fmap selected)
  . apollonianTrees
  $ cs
  where
    p c = radius c >= thresh

-- | Given a set of /four/ mutually tangent circles, generate the
--   infinite Apollonian tree rooted at the given set, represented as
--   a list of four subtrees.  Each node in the tree is a kissing set
--   with one circle selected which has just been flipped.  The three
--   children of a node represent the kissing sets obtained by
--   selecting each of the other three circles and flipping them.  The
--   initial roots of the four trees are chosen by selecting and
--   flipping each of the circles in the starting set. This
--   representation has the property that each circle in the
--   Apollonian gasket is the selected circle in exactly one node
--   (except that the initial four circles never appear as the
--   selected circle in any node).
apollonianTrees :: RealFloat n => [Circle n] -> [Tree (KissingSet (Circle n))]
apollonianTrees = map (apollonianTree . flipSelected) . kissingSets

-- | Generate a single Apollonian tree from a root kissing set.  See
--   the documentation for 'apollonianTrees' for an explanation.
apollonianTree :: RealFloat n => KissingSet (Circle n) -> Tree (KissingSet (Circle n))
apollonianTree = unfoldTree (id &&& (map flipSelected . selectOthers))

-- | Prune a tree at the shallowest points where the predicate is not
--   satisfied.
prune :: (a -> Bool) -> Tree a -> Maybe (Tree a)
prune p (Node a ts)
  | not (p a) = Nothing
  | otherwise = Just $ Node a (catMaybes (map (prune p) ts))

------------------------------------------------------------
--  Diagram generation
------------------------------------------------------------

-- | Draw a circle.
drawCircle :: (Renderable (Path V2 n) b, TypeableFloat n) =>
              Circle n -> QDiagram b V2 n Any
drawCircle c = circle (radius c) # moveTo (center c)
                                 # fcA transparent

-- | Draw a generated gasket, using a line width 0.003 times the
--   radius of the largest circle.
drawGasket :: (Renderable (Path V2 n) b, TypeableFloat n) =>
              [Circle n] -> QDiagram b V2 n Any
drawGasket cs = F.foldMap drawCircle cs

-- | Draw an Apollonian gasket: the first argument is the threshold;
--   the recursion will stop upon reaching circles with radii less than
--   it. The next three arguments are bends of three circles.
apollonianGasket :: (Renderable (Path V2 n) b, TypeableFloat n)
                 => n -> n -> n -> n -> QDiagram b V2 n Any
apollonianGasket thresh b1 b2 b3 = drawGasket . apollonian thresh $ (initialConfig b1 b2 b3)

------------------------------------------------------------
-- Some notes on floating-point error
-- (only for the intrepid)
------------------------------------------------------------
{-

-- code from Gerald Gutierrez, personal communication

module Main where

import           Data.Complex
import           Diagrams.Backend.SVG.CmdLine
import           Diagrams.Prelude

-- ------^---------^---------^---------^---------^---------^---------^--------

data Circle = Circle Double (Complex Double) deriving (Show)

descartes a b c
  = (s + r, s - r)
  where
    s = a + b + c
    r = 2 * sqrt ((a * b) + (b * c) + (c * a))

descartesDual a b c d
  = 2 * (a + b + c) - d

soddies (Circle k1 b1) (Circle k2 b2) (Circle k3 b3)
  = ( Circle k4 b4
    , Circle k5 b5 )
  where
    (k4, k5) = descartes k1 k2 k3
    (b4, b5) = descartes b1 b2 b3

soddiesDual (Circle k1 b1) (Circle k2 b2) (Circle k3 b3) (Circle k4 b4)
  = Circle (descartesDual k1 k2 k3 k4) (descartesDual b1 b2 b3 b4)

mutuallyTangentCirclesFromTriangle z1 z2 z3
  = ( Circle k1 (z1 * (k1 :+ 0))
    , Circle k2 (z2 * (k2 :+ 0))
    , Circle k3 (z3 * (k3 :+ 0)) )
  where
    a = magnitude (z2 - z3)
    b = magnitude (z3 - z1)
    c = magnitude (z1 - z2)
    s = (a + b + c) / 2
    k1 = 1 / (s - a)
    k2 = 1 / (s - b)
    k3 = 1 / (s - c)

main :: IO ()
main = mainWith picture

pic = mainWith picture

mkCircle :: Circle -> Diagram B
mkCircle (Circle k b)
  = circle (1 / k) # moveTo (p2 (realPart z, imagPart z))
  where
    z = b / (k :+ 0)

picture :: Diagram B
picture
  = mkCircle c1 <> mkCircle c2 <> mkCircle c3 <> mkCircle c4 <> mkCircle c5
  where
    (c1, c2, c3) = mutuallyTangentCirclesFromTriangle z1 z2 z3
    (c4, c5) = soddies c1 c2 c3

    -- z1 = 0 :+ 0
    -- z2 = 3 :+ 0
    -- z3 = 0 :+ 4

    -- z1 = (-0.546) :+ (-0.755)
    -- z2 = ( 0.341) :+ (-0.755)
    -- z3 = (-0.250) :+ ( 0.428)

ofsBad = 0.15397 -- doesn't work
ofsGood = 0.15398 -- works

ofs = ofsGood

z1 = ofs + ((-0.546) :+ (-0.755))
z2 = ofs + (( 0.341) :+ (-0.755))
z3 = ofs + ((-0.250) :+ ( 0.428))

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

Email to Gerald Gutierrez, 30 Sep 2016:

I got a chance to sit down and think hard about your code today, and
I think I have figured out what the problem is.

If you look at the outputs of 'soddies c1 c2 c3' using the two
different values for 'ofs', you will see that they are *almost* the
same, except that the complex numbers have been switched (though in
fact their real components are almost the same so you only notice the
difference in the imaginary components).  This clearly causes
incorrect results since the y-coordinates represented by the imaginary
components are now getting scaled by different bend values.  So the
resulting circles are of the right size and have (close to) the
correct x-coordinates, but their y-coordinates are wrong.

The problem is that in the 'soddies' function, you make independent
calls to 'descartes k1 k2 k3' (to solve for the bends) and 'descartes
b1 b2 b3' (to solve for the bend-center products), and then *assume*
that they return their solutions *in the same order*, so that you can
match up the first values to make one circle and the second values to
make another circle.  I would guess that this is *usually* true but
apparently not when certain values are hovering right around a branch
cut of sqrt, or something like that.

I think my code in Diagrams.TwoD.Apollonian suffers from the same
problem.  Ultimately I think it is due to the inherent inaccuracy of
floating-point numbers; there is nothing wrong with the formula
itself.  It would be nice to figure out how to correct for this, but I
am not sure how off the top of my head.  The interesting thing is that
switching the bend-center products does not seem to violate the
generalized Descartes' Theorem at all --- so it would seem that it
does not actually completely characterize mutually tangent circles,
that is, there exist sets of circles satisfying the theorem which are
not in fact mutually tangent.

-}

max_stars_repo_name: phadej/diagrams-contrib

max_issues_repo_name: phadej/diagrams-contrib

max_forks_repo_name: phadej/diagrams-contrib


text:
{-# LANGUAGE BangPatterns #-}
-- # LANGUAGE UndecidableInstances #
{-# LANGUAGE GADTs #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
module Data.Network.FullyConnected
  (FeedForwardFC, FullyConnected(FC), Layers,
   weights, biases, activation,
   unGrad, unGradBatch, Grad(GradBatch, Grad))
where

import Numeric.LinearAlgebra as LA hiding (fromList, toList)
import Numeric.LinearAlgebra.Devel
import qualified Data.Vector.Storable as V hiding (foldr, modify)
import System.Random
import Control.DeepSeq
-- import Data.Foldable as F
import Data.Bifunctor
import Data.Network as Network
import Data.Network.Activation
import Algorithm.Training
import Algorithm.GradientDescent
import Data.Network.Utils (squeeze, (.*))

-- Fully connected feedforward layer
data FullyConnected i a = FC {weights :: {-# UNPACK #-} !(Matrix a),
                              biases ::  {-# UNPACK #-} !(Vector a),
                              activation :: ActivationFunction i
                       }

instance (V.Storable a, NFData a) => NFData (FullyConnected i a) where
  rnf (FC w b a) = rnf w `deepseq` rnf b

instance Layer (Vector R) (FullyConnected (Vector R) R) where
  compute (FC weights biases _) input = biases + (input <# weights)
  run     l@(FC weights biases af) input = unAct af (compute l input)
  makeRandomLayer g rf i o af = let rnds = (randomRs (0,1) g) :: [Double]
                                    rnds' = (rf rnds)
                                    rw   = (i><o) rnds'
                                    rb   = o |> (repeat 0)
                                in (FC rw rb af)

instance Layer (Matrix R) (FullyConnected (Matrix R) R) where
  compute (FC weights biases _) input = (fromRows $ replicate (rows input) biases) + (input <> weights)
  run     l@(FC _ _ af) input = unAct af (compute l input)
  makeRandomLayer g rf i o af = let rnds = (randomRs (0,1) g) :: [Double]
                                    rnds' = (rf rnds)
                                    rw   = (i><o) rnds'
                                    rb   = o |> (repeat 0)
                                in (FC rw rb af)

instance (Container Matrix a, Container Vector a, Num a) => Show (FullyConnected i a) where
  show l = let (a,b) = size $ weights l
           in "Layer " ++ (show a) ++ "x" ++ (show b) ++ " " ++ (name $ activation l)

instance Floating a => Num (FullyConnected a Double) where
  (FC w b a) + (FC w' b' _) = FC (w + w') (b + b') a
  (FC w b a) - (FC w' b' _) = FC (w - w') (b - b') a
  (FC w b a) * (FC w' b' _) = FC (w * w') (b * b') a
  abs (FC w b a) = FC (abs w) (abs b) a
  signum (FC w b a) = FC (signum w) (signum b) a
  fromInteger x = FC (fromInteger x) (fromInteger x) (act id "identity")

instance Floating a => Fractional (FullyConnected a Double) where
  (FC w b a) / (FC w' b' _) = FC (w / w') (b / b') a
  recip (FC w b a) = FC (recip w) (recip b) a
  fromRational a = FC (fromRational a) (fromRational a) (act id "identity")

instance (Container Vector a, Container Matrix a) => Linear a (FullyConnected i) where
  scale x (FC w b a) = FC (scale x w) (scale x b) a

-- Container Layer instance?

instance (Layer i (FullyConnected i Double)) => Network i (Layers (FullyConnected i Double)) where
  runNetwork (Output l) i = run l i
  runNetwork (l :=> n)  i = runNetwork n (fst $ run l i)
  runLayers  (Output l) i = [run l i]
  runLayers  (l :=> n)  i = let y = (run l i)
                            in  y:(runLayers n (fst y))

  makeRandomNetwork g t nt i [] (o, af) = let l = makeRandomLayer g t i o af
                                          in  (Output l)
  makeRandomNetwork g t nt i ((h,haf):hs) o = let (g1,g2) = split g
                                                  l = makeRandomLayer g1 t i h haf
                                              in nt $ l :=> (makeRandomNetwork g2 t nt h hs o)

type FeedForwardFC i = Layers (FullyConnected i Double)

instance (Monad m) => Backpropagation Data (NetworkEvaluator Data) (FeedForwardFC Data) m where
  backward lf n (out,tar) das = do
    let δout = (derivate lf (tar, out)) -- dE/dy
        deltas = scanr (\(l, a') δ -> let w = weights l in a' * (w #> δ)) δout (zip (tail $ toList n) das)
    return (deltas) -- deltas for computing bias

instance (Monad m) => Backpropagation Batch (NetworkEvaluator Batch) (FeedForwardFC Batch) m where
  backward lf n (out,tar) das = do
    let δout = tr (derivate lf (tar, out)) -- dE/dy
        deltas = scanr (\(l, a') δ -> let w = weights l in (tr a') * (w <> δ)) δout (zip (tail $ toList n) das)
    return (deltas) -- deltas for computing bias

-- type instance Params lf = forall i. ((FeedForwardFC i), (i,i))

instance (Monad m) => GradientDescent (NetworkEvaluator Batch) ((FeedForwardFC Batch), (Batch,Batch)) m where
  -- parameters of loss function: Network(weights, biases, activation functions for each layer)
  -- and sample((input, output) pair)

  data Grad (NetworkEvaluator Batch) = GradBatch { unGradBatch :: ([Matrix Double], [Vector Double]) } deriving (Show)

  -- Gradients over parameters(weights, biases)
  grad lf (n, (i,t)) = do
    let (as, as') = unzip $ runLayers n i -- forward propagation: compute layers outputs and activation derivatives
        (out) = last as
    (ds) <- backward lf n (out, t) (init as') -- compute deltas with backpropagation
    let r  = fromIntegral $ rows i -- size of minibatch
    let gs = zipWith (\δ a -> tr (δ <> a)) ds (i:init as)
    return $ GradBatch ((recip r .*) <$> gs, ((recip r .*) . squeeze) <$> ds)

  move lr (n, (i,t)) (GradBatch (gs, ds)) = do
    -- update function
    -- divide by batch size
    let update = (\(FC w b af) g δ -> FC (w + (lr).*g) (b + (lr).*δ) af)
        n' = Network.fromList $ zipWith3 update (Network.toList n) gs ds -- new model
    return (n', (i,t))

instance (Monad m) => GradientDescent (NetworkEvaluator Data) ((FeedForwardFC Data), (Data,Data)) m where
  -- parameters of loss function: Network(weights, biases, activation functions for each layer)
  -- and sample((input, output) pair)

  data Grad (NetworkEvaluator Data) = Grad { unGrad :: ([Matrix Double], [Vector Double]) } deriving (Show)

  -- Gradients over parameters(weights, biases)
  grad lf (n, (i,t)) = do
    let aos = runLayers n i -- forward propagation: compute layers outputs and activation derivatives
        (as, as') = (unzip) aos -- juggle data to extract layers inputs and derivatives in right order
        out = last as
    (ds) <- backward lf n (out,t) (init as') -- compute gradients with backpropagation
    let gs = zipWith (\δ i -> tr (δ `outer` i)) ds (i:as)
    return $ Grad (gs, ds)

  move lr (n, (i,t)) (Grad (gs, ds)) = do
    -- update function
    let update = (\(FC w b af) g δ -> FC (w + lr.*g) (b + lr.* δ) af)
        n' = Network.fromList $ zipWith3 update (Network.toList n) gs ds -- new model
    return (n', (i,t))

max_stars_repo_name: DrPyser/NeuralNetwork

max_issues_repo_name: DrPyser/NeuralNetwork

max_forks_repo_name: DrPyser/NeuralNetwork


text:
module AI.Learning.Core where

import Numeric.LinearAlgebra
import Numeric.GSL.Minimization

---------------
-- Utilities --
---------------

-- |Sigmoid function:
--
--  > sigmoid x = 1 / (1 + exp (-x))
--
--  Used in the logistic regression and neural network modules.
sigmoid :: Floating a => a -> a
sigmoid x = 1 / (1 + exp (-x))

------------------
-- Optimization --
------------------

-- |Simplified minimization function. You supply functions that compute the
--  quantity to be minimized and its gradient, and an initial guess, and the
--  final solution is returned.
minimizeS :: (Vector Double -> Double)          -- f
          -> (Vector Double -> Vector Double)   -- gradient
          -> Vector Double                      -- initial x
          -> Vector Double
minimizeS f g x = fst $ minimizeVD VectorBFGS2 prec niter sz tol f g x
    where prec    = 1e-9
          niter   = 1000
          sz      = 0.1
          tol     = 0.1

max_stars_repo_name: cagix/aima-haskell


text:
{-# LANGUAGE AllowAmbiguousTypes, CPP, ConstraintKinds, DataKinds,
             FlexibleInstances, GADTs, MultiParamTypeClasses, PolyKinds,
             QuantifiedConstraints, TypeFamilyDependencies, TypeOperators,
             UndecidableInstances #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}
-- | The class defining dual components of dual numbers and related classes,
-- type families, constraints and instances. This is a low-level API
-- used to define types and operations in "HordeAd.Core.DualNumber"
-- that is the high-level API.
module HordeAd.Core.DualClass
  ( IsPrimalWithScalar, IsPrimalAndHasFeatures, IsScalar, HasDelta
  , DMode(..), Dual, IsPrimal(..), HasRanks(..)
  , HasVariables(..)  -- use sparringly
  ) where

import Prelude

import qualified Data.Array.Convert
import qualified Data.Array.Dynamic as OTB
import qualified Data.Array.DynamicS as OT
import qualified Data.Array.Shaped as OSB
import qualified Data.Array.ShapedS as OS
import           Data.MonoTraversable (Element, MonoFunctor)
import           Data.Proxy (Proxy)
import qualified Data.Strict.Vector as Data.Vector
import qualified Data.Vector.Generic as V
import           GHC.TypeLits (KnownNat, natVal, type (+))
import           Numeric.LinearAlgebra (Matrix, Numeric, Vector)
import qualified Numeric.LinearAlgebra as HM

import HordeAd.Internal.Delta

-- * Abbreviations to export (not used anywhere below)

-- | The intended semantics (not fully enforced by the constraint in isolation)
-- is that the second type is the primal component of a dual number type
-- at an unknown rank, with the given differentiation mode
-- and underlying scalar.
type IsPrimalWithScalar (d :: DMode) a r =
  (ScalarOf a ~ r, IsPrimal d a, HasVariables a)

-- | A shorthand for a useful set of constraints.
type IsPrimalAndHasFeatures (d :: DMode) a r =
  (IsPrimalWithScalar d a r, RealFloat a, MonoFunctor a, Element a ~ r)

-- | A mega-shorthand for a bundle of connected type constraints.
-- The @Scalar@ in the name means that the second argument is the underlying
-- scalar type of a well behaved (wrt the differentiation mode in the first
-- argument) collection of primal and dual components of dual numbers.
type IsScalar (d :: DMode) r =
  ( HasRanks d r, Ord r, Numeric r, Show r
  , IsPrimalAndHasFeatures d r r
  , IsPrimalAndHasFeatures d (Vector r) r
  , IsPrimalAndHasFeatures d (Matrix r) r
  , IsPrimalAndHasFeatures d (OT.Array r) r
  -- This fragment is for @OS.Array@ and it's irregular, because we can't
  -- mention @sh@ and so fully apply the type constructor.
  , IsPrimalS d r  -- TODO: Floating (OS.Array sh r), MonoFunctor
  )

-- | Is a scalar and will be used to compute gradients via delta-expressions.
type HasDelta r = ( IsScalar 'DModeGradient r
                  , Dual 'DModeGradient r ~ Delta0 r )

-- * Class definitions

-- | The enumeration of all possible differentiation (and more generally,
-- computation with dual numbers) schemes.
data DMode =
    DModeGradient
  | DModeDerivative

-- | The type family that enumerates all possible "ranks"
-- for each differentiation mode.
-- The second type argument is meant to be the primal component
-- of dual numbers. The result is the dual component.
--
-- Rank 0 is special because, in derivatives mode, the dual component
-- is not the primal component wrapped in a datatype or newtype constructor.
-- This makes impossible a representation of primal and dual components as
-- the primal plus the type constructor for creating the dual.
--
-- Rank S is special, because of the extra type parameter @sh@ representing
-- a shape. This is another obstacle to a dual number representation via
-- a single-argument type constructor.
type family Dual (d :: DMode) a = result | result -> d a where
  Dual 'DModeGradient Double = Delta0 Double
  Dual 'DModeGradient Float = Delta0 Float
  Dual 'DModeGradient (Vector r) = Delta1 r
  Dual 'DModeGradient (Matrix r) = Delta2 r
  Dual 'DModeGradient (OT.Array r) = DeltaX r
  Dual 'DModeGradient (OS.Array sh r) = DeltaS sh r
-- not injective:  Dual 'DModeDerivative r = r
  Dual 'DModeDerivative Double = Double
  Dual 'DModeDerivative Float = Float
  Dual 'DModeDerivative (Vector r) = Vector r
  Dual 'DModeDerivative (Matrix r) = Matrix r
  Dual 'DModeDerivative (OT.Array r) = OT.Array r
  Dual 'DModeDerivative (OS.Array sh r) = OS.Array sh r

-- | The underlying scalar of a given primal component of a dual number.
-- A long name to remember not to use, unless necessary, and not to export.
type family ScalarOf a where
  ScalarOf Double = Double
  ScalarOf Float = Float
  ScalarOf (Vector r) = r
  ScalarOf (Matrix r) = r
  ScalarOf (OT.Array r) = r
  ScalarOf (OS.Array sh r) = r

-- | Second argument is a primal component of dual numbers at some rank
-- wrt the differentiation mode given in the first argument.
class IsPrimal d a where
  dZero :: Dual d a
  dScale :: a -> Dual d a -> Dual d a
  dAdd :: Dual d a -> Dual d a -> Dual d a
  dDelay :: Dual d a -> Dual d a

-- | Part 1/2 of a hack to squeeze the shaped tensors rank,
-- with its extra @sh@ parameter, into the 'IsPrimal' class.
class IsPrimalS d r where
  dZeroS :: forall sh. OS.Shape sh => Dual d (OS.Array sh r)
  dScaleS :: forall sh. OS.Shape sh
          => OS.Array sh r -> Dual d (OS.Array sh r) -> Dual d (OS.Array sh r)
  dAddS :: forall sh. OS.Shape sh
        => Dual d (OS.Array sh r) -> Dual d (OS.Array sh r)
        -> Dual d (OS.Array sh r)
  dDelayS :: forall sh. OS.Shape sh
             => Dual d (OS.Array sh r) -> Dual d (OS.Array sh r)

-- | Part 2/2 of a hack to squeeze the shaped tensors rank,
-- with its extra @sh@ parameter, into the 'IsPrimal' class.
instance (IsPrimalS d r, OS.Shape sh) => IsPrimal d (OS.Array sh r) where
  dZero = dZeroS
  dScale = dScaleS
  dAdd = dAddS
  dDelay = dDelayS

-- | Assuming that the first argument is the primal component of dual numbers
-- with the underyling scalar in the second argument and with differentiation
-- mode `DModeGradient`, it additionally admits delta-variable
-- introduction and binding as defined by the methods of the class.
class HasVariables a where
  dVar :: DeltaId a -> Dual 'DModeGradient a
  bindInState :: Dual 'DModeGradient a
              -> DeltaState (ScalarOf a)
              -> (DeltaState (ScalarOf a), DeltaId a )
  dOutline :: CodeOut -> [a] -> [Dual 'DModeGradient a]
           -> Dual 'DModeGradient a

-- | The class provides methods required for the second type parameter
-- to be the underlying scalar of a well behaved collection of dual numbers
-- of various ranks wrt the differentation mode given in the first argument.
class HasRanks (d :: DMode) r where
  dSumElements0 :: Dual d (Vector r) -> Int -> Dual d r
  dIndex0 :: Dual d (Vector r) -> Int -> Int -> Dual d r
  dDot0 :: Vector r -> Dual d (Vector r) -> Dual d r
  dFromX0 :: Dual d (OT.Array r) -> Dual d r
  dFromS0 :: Dual d (OS.Array '[] r) -> Dual d r

  dSeq1 :: Data.Vector.Vector (Dual d r) -> Dual d (Vector r)
  dKonst1 :: Dual d r -> Int -> Dual d (Vector r)
  dAppend1 :: Dual d (Vector r) -> Int -> Dual d (Vector r) -> Dual d (Vector r)
  dSlice1 :: Int -> Int -> Dual d (Vector r) -> Int -> Dual d (Vector r)
  dSumRows1 :: Dual d (Matrix r) -> Int -> Dual d (Vector r)
  dSumColumns1 :: Dual d (Matrix r) -> Int -> Dual d (Vector r)
  dM_VD1 :: Matrix r -> Dual d (Vector r) -> Dual d (Vector r)
  dMD_V1 :: Dual d (Matrix r) -> Vector r -> Dual d (Vector r)
  dFromX1 :: Dual d (OT.Array r) -> Dual d (Vector r)
  dFromS1 :: KnownNat len
          => Dual d (OS.Array '[len] r) -> Dual d (Vector r)
  dReverse1 :: Dual d (Vector r) -> Dual d (Vector r)
  dFlatten1 :: Int -> Int -> Dual d (Matrix r) -> Dual d (Vector r)
  dFlattenX1 :: OT.ShapeL -> Dual d (OT.Array r) -> Dual d (Vector r)
  dFlattenS1 :: OS.Shape sh
             => Dual d (OS.Array sh r) -> Dual d (Vector r)

  dFromRows2 :: Data.Vector.Vector (Dual d (Vector r)) -> Dual d (Matrix r)
  dFromColumns2 :: Data.Vector.Vector (Dual d (Vector r)) -> Dual d (Matrix r)
  dKonst2 :: Dual d r -> (Int, Int) -> Dual d (Matrix r)
  dTranspose2 :: Dual d (Matrix r) -> Dual d (Matrix r)
  dM_MD2 :: Matrix r -> Dual d (Matrix r) -> Dual d (Matrix r)
  dMD_M2 :: Dual d (Matrix r) -> Matrix r -> Dual d (Matrix r)
  dRowAppend2 :: Dual d (Matrix r) -> Int -> Dual d (Matrix r) -> Dual d (Matrix r)
  dColumnAppend2 :: Dual d (Matrix r) -> Int -> Dual d (Matrix r) -> Dual d (Matrix r)
  dRowSlice2 :: Int -> Int -> Dual d (Matrix r) -> Int -> Dual d (Matrix r)
  dColumnSlice2 :: Int -> Int -> Dual d (Matrix r) -> Int -> Dual d (Matrix r)
  dAsRow2 :: Dual d (Vector r) -> Dual d (Matrix r)
  dAsColumn2 :: Dual d (Vector r) -> Dual d (Matrix r)
  dFromX2 :: Dual d (OT.Array r) -> Dual d (Matrix r)
  dFromS2 :: (KnownNat rows, KnownNat cols)
          => Dual d (OS.Array '[rows, cols] r) -> Dual d (Matrix r)

  dFlipud2 :: Dual d (Matrix r) -> Dual d (Matrix r)
  dFliprl2 :: Dual d (Matrix r) -> Dual d (Matrix r)
  dReshape2 :: Int -> Dual d (Vector r) -> Dual d (Matrix r)
  dConv2 :: Matrix r -> Dual d (Matrix r) -> Dual d (Matrix r)

  dKonstX :: Dual d r -> OT.ShapeL -> Dual d (OT.Array r)
  dAppendX :: Dual d (OT.Array r) -> Int -> Dual d (OT.Array r) -> Dual d (OT.Array r)
  dSliceX :: Int -> Int -> Dual d (OT.Array r) -> Int -> Dual d (OT.Array r)
  dIndexX :: Dual d (OT.Array r) -> Int -> Int -> Dual d (OT.Array r)
  dRavelFromListX :: [Dual d (OT.Array r)] -> Dual d (OT.Array r)
  dReshapeX :: OT.ShapeL -> OT.ShapeL -> Dual d (OT.Array r) -> Dual d (OT.Array r)
  dFrom0X :: Dual d r -> Dual d (OT.Array r)
  dFrom1X :: Dual d (Vector r) -> Dual d (OT.Array r)
  dFrom2X :: Dual d (Matrix r) -> Int -> Dual d (OT.Array r)
  dFromSX :: OS.Shape sh
          => Dual d (OS.Array sh r) -> Dual d (OT.Array r)

  dKonstS :: OS.Shape sh
          => Dual d r -> Dual d (OS.Array sh r)
  dAppendS :: (OS.Shape sh, KnownNat m, KnownNat n)
           => Dual d (OS.Array (m ': sh) r) -> Dual d (OS.Array (n ': sh) r)
           -> Dual d (OS.Array ((m + n) ': sh) r)
  dSliceS :: (KnownNat i, KnownNat n, KnownNat k, OS.Shape rest)
          => Proxy i -> Proxy n -> Dual d (OS.Array (i + n + k ': rest) r)
          -> Dual d (OS.Array (n ': rest) r)
  dIndexS :: (KnownNat ix, KnownNat k, OS.Shape rest)
          => Dual d (OS.Array (ix + 1 + k ': rest) r)
          -> Proxy ix -> Dual d (OS.Array rest r)
  dRavelFromListS :: (KnownNat k, OS.Shape rest)
                  => [Dual d (OS.Array rest r)]
                  -> Dual d (OS.Array (k : rest) r)
  dReshapeS :: (OS.Shape sh, OS.Shape sh', OS.Size sh ~ OS.Size sh')
            => Dual d (OS.Array sh r) -> Dual d (OS.Array sh' r)
  dFrom0S :: Dual d r -> Dual d (OS.Array '[] r)
  dFrom1S :: KnownNat n => Dual d (Vector r) -> Dual d (OS.Array '[n] r)
  dFrom2S :: (KnownNat rows, KnownNat cols)
          => Proxy cols
          -> Dual d (Matrix r) -> Dual d (OS.Array '[rows, cols] r)
  dFromXS :: OS.Shape sh => Dual d (OT.Array r) -> Dual d (OS.Array sh r)

-- * Backprop gradient method instances

instance IsPrimal 'DModeGradient Double where
  dZero = Zero0
  dScale = Scale0
  dAdd = Add0
  dDelay = Delay0

instance IsPrimal 'DModeGradient Float where
  -- Identical as above:
  dZero = Zero0
  dScale = Scale0
  dAdd = Add0
  dDelay = Delay0

instance IsPrimal 'DModeGradient (Vector r) where
  dZero = Zero1
  dScale = Scale1
  dAdd = Add1
  dDelay = Delay1

instance IsPrimal 'DModeGradient (Matrix r) where
  dZero = Zero2
  dScale = Scale2
  dAdd = Add2
  dDelay = Delay2

instance IsPrimal 'DModeGradient (OT.Array r) where
  dZero = ZeroX
  dScale = ScaleX
  dAdd = AddX
  dDelay = DelayX

instance IsPrimalS 'DModeGradient r where
  dZeroS = ZeroS
  dScaleS = ScaleS
  dAddS = AddS
  dDelayS = DelayS

instance HasVariables Double where
  dVar = Var0
  {-# INLINE bindInState #-}
  bindInState = bindInState0
  dOutline = Outline0

instance HasVariables Float where
  dVar = Var0
  {-# INLINE bindInState #-}
  bindInState = bindInState0
  dOutline = Outline0

instance HasVariables (Vector r) where
  dVar = Var1
  {-# INLINE bindInState #-}
  bindInState = bindInState1
  dOutline = Outline1

instance HasVariables (Matrix r) where
  dVar = Var2
  {-# INLINE bindInState #-}
  bindInState = bindInState2
  dOutline = Outline2

instance HasVariables (OT.Array r) where
  dVar = VarX
  {-# INLINE bindInState #-}
  bindInState = bindInStateX
  dOutline = OutlineX

instance OS.Shape sh => HasVariables (OS.Array sh r) where
  dVar = VarS
  {-# INLINE bindInState #-}
  bindInState u' st = let (st2, did) = bindInStateX (FromSX u') st
                      in (st2, convertDeltaId did)
  dOutline = OutlineS

instance Dual 'DModeGradient r ~ Delta0 r
         => HasRanks 'DModeGradient r where
  dSumElements0 = SumElements0
  dIndex0 = Index0
  dDot0 = Dot0
  dFromX0 = FromX0
  dFromS0 = FromS0
  dSeq1 = Seq1
  dKonst1 = Konst1
  dAppend1 = Append1
  dSlice1 = Slice1
  dSumRows1 = SumRows1
  dSumColumns1 = SumColumns1
  dM_VD1 = M_VD1
  dMD_V1 = MD_V1
  dFromX1 = FromX1
  dFromS1 = FromS1
  dReverse1 = Reverse1
  dFlatten1 = Flatten1
  dFlattenX1 = FlattenX1
  dFlattenS1 = FlattenS1
  dFromRows2 = FromRows2
  dFromColumns2 = FromColumns2
  dKonst2 = Konst2
  dTranspose2 = Transpose2
  dM_MD2 = M_MD2
  dMD_M2 = MD_M2
  dRowAppend2 = RowAppend2
  dColumnAppend2 = ColumnAppend2
  dRowSlice2 = RowSlice2
  dColumnSlice2 = ColumnSlice2
  dAsRow2 = AsRow2
  dAsColumn2 = AsColumn2
  dFromX2 = FromX2
  dFromS2 = FromS2
  dFlipud2 = Flipud2
  dFliprl2 = Fliprl2
  dReshape2 = Reshape2
  dConv2 = Conv2
  dKonstX = KonstX
  dAppendX = AppendX
  dSliceX = SliceX
  dIndexX = IndexX
  dRavelFromListX = RavelFromListX
  dReshapeX = ReshapeX
  dFrom0X = From0X
  dFrom1X = From1X
  dFrom2X = From2X
  dFromSX = FromSX
  dKonstS = KonstS
  dAppendS = AppendS
  dSliceS = SliceS
  dIndexS = IndexS
  dRavelFromListS = RavelFromListS
  dReshapeS = ReshapeS
  dFrom0S = From0S
  dFrom1S = From1S
  dFrom2S = From2S
  dFromXS = FromXS

-- * Alternative instances: forward derivatives computed on the spot

instance IsPrimal 'DModeDerivative Double where
  dZero = 0
  dScale k d = k * d
  dAdd d e = d + e
  dDelay = id  -- no delaying

instance IsPrimal 'DModeDerivative Float where
  dZero = 0
  dScale k d = k * d
  dAdd d e = d + e
  dDelay = id

-- These constraints force @UndecidableInstances@.
instance Num (Vector r)
         => IsPrimal 'DModeDerivative (Vector r) where
  dZero = 0
  dScale k d = k * d
  dAdd d e = d + e
  dDelay = id

instance Num (Matrix r)
         => IsPrimal 'DModeDerivative (Matrix r) where
  dZero = 0
  dScale k d = k * d
  dAdd d e = d + e
  dDelay = id

instance Num (OT.Array r)
         => IsPrimal 'DModeDerivative (OT.Array r) where
  dZero = 0
  dScale k d = k * d
  dAdd d e = d + e
  dDelay = id

instance (Numeric r, Num (Vector r))
         => IsPrimalS 'DModeDerivative r where
  dZeroS = 0
  dScaleS k d = k * d
  dAddS d e = d + e
  dDelayS = id

instance ( Numeric r, Num (Vector r)
         , Dual 'DModeDerivative r ~ r )
         => HasRanks 'DModeDerivative r where
  dSumElements0 vd _ = HM.sumElements vd
  dIndex0 d ix _ = d V.! ix
  dDot0 = (HM.<.>)
  dFromX0 = OT.unScalar
  dFromS0 = OS.unScalar
  dSeq1 = V.convert
  dKonst1 = HM.konst
  dAppend1 d _k e = d V.++ e
  dSlice1 i n d _len = V.slice i n d
  dM_VD1 = (HM.#>)
  dMD_V1 = (HM.#>)
  dSumRows1 dm _cols = V.fromList $ map HM.sumElements $ HM.toRows dm
  dSumColumns1 dm _rows = V.fromList $ map HM.sumElements $ HM.toColumns dm
  dFromX1 = OT.toVector
  dFromS1 = OS.toVector
  dReverse1 = V.reverse
  dFlatten1 _rows _cols = HM.flatten
  dFlattenX1 _sh = OT.toVector
  dFlattenS1 = OS.toVector
  dFromRows2 = HM.fromRows . V.toList
  dFromColumns2 = HM.fromColumns . V.toList
  dKonst2 = HM.konst
  dTranspose2 = HM.tr'
  dM_MD2 = (HM.<>)
  dMD_M2 = (HM.<>)
  dAsRow2 = HM.asRow
  dAsColumn2 = HM.asColumn
  dRowAppend2 d _k e = d HM.=== e
  dColumnAppend2 d _k e = d HM.||| e
  dRowSlice2 i n d _rows = HM.takeRows n $ HM.dropRows i d
  dColumnSlice2 i n d _cols = HM.takeColumns n $ HM.dropColumns i d
  dFromX2 d = case OT.shapeL d of
    [_rows, cols] -> HM.reshape cols $ OT.toVector d
    _ -> error "dFromX2: wrong tensor dimensions"
  dFromS2 d = case OS.shapeL d of
    [_rows, cols] -> HM.reshape cols $ OS.toVector d
    _ -> error "dFromS2: wrong tensor dimensions"
  dFlipud2 = HM.flipud
  dFliprl2 = HM.fliprl
  dReshape2 = HM.reshape
  dConv2 = HM.conv2
  dKonstX d sz = OT.constant sz d
  dAppendX d _k e = d `OT.append` e
  dSliceX i n d _len = OT.slice [(i, n)] d
  dIndexX d ix _len = OT.index d ix
  dRavelFromListX ld =
    let sh = case ld of
          d : _ -> length ld : OT.shapeL d
          [] -> []
    in OT.ravel $ OTB.fromList sh ld
  dReshapeX _sh = OT.reshape
  dFrom0X = OT.scalar
  dFrom1X d = OT.fromVector [V.length d] d
  dFrom2X d cols = OT.fromVector [HM.rows d, cols] $ HM.flatten d
  dFromSX = Data.Array.Convert.convert
  dKonstS = OS.constant
  dAppendS = OS.append
  dSliceS (_ :: Proxy i) (_ :: Proxy n) = OS.slice @'[ '(i, n) ]
  dIndexS d proxyIx = OS.index d (fromInteger $ natVal proxyIx)
  dRavelFromListS = OS.ravel . OSB.fromList
  dReshapeS = OS.reshape
  dFrom0S = OS.scalar
  dFrom1S = OS.fromVector
  dFrom2S _ = OS.fromVector . HM.flatten
  dFromXS = Data.Array.Convert.convert


text:
{-# LANGUAGE BangPatterns #-}

module Main where
--import qualified Data.Vector.Storable as V
import Foreign
import Data.Ratio
import Data.Complex

complexZI :: Complex Int
complexZI = 1 :+ 1

ratio23 :: Ratio Int
ratio23 =  1 % 1

putter :: Storable a => a -> Ptr a -> IO a
putter v !ptr = do poke ptr v ; peek ptr

main =
  do
      !vComplex <- alloca (putter complexZI)
      !vRatio <- alloca (putter ratio23)
      if vComplex  == complexZI &&  vRatio  == ratio23
        then putStrLn "success"
        else putStrLn "uh oh, something is wrong with storable"

max_stars_repo_name: JakobBruenker/ghc

max_issues_repo_name: JakobBruenker/ghc

max_forks_repo_name: JakobBruenker/ghc