text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-}

{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedLabels      #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-|
Module      : Grenade.Layers.BatchNormalisation
Description : Batch normalization layer
Copyright   : (c) Theo Charalambous, 2020
License     : BSD2
Stability   : experimental

The layer follows the implementation as described in the paper:
Sergey Ioffe and Christian Szegedy. Batch normalization:  Accelerating deep network training by reducing internal covariate shift
<http://arxiv.org/abs/1502.03167>
-}

module Grenade.Layers.BatchNormalisation
  (
  -- * Layer Definition
    BatchNorm (..)
  , BatchNormParams(..)

  -- * Differentation Types
  , BatchNormTape(..)
  , BatchNormGrad(..)

  -- * Helper functions
  , initBatchNorm
  )
where

import           Control.DeepSeq
import           Data.Kind                         (Type)
import           Data.List                         (transpose, zipWith5)
import           Data.Maybe                        (fromJust)
import           Data.Proxy
import           Data.Serialize
import           GHC.TypeLits

import qualified Numeric.LinearAlgebra             as LA
import           Numeric.LinearAlgebra.Static      hiding (Seed)

import           Grenade.Core
import           Grenade.Layers.Internal.BatchNorm
import           Grenade.Layers.Internal.Update
import           Grenade.Onnx
import           Grenade.Types
import           Grenade.Utils.LinearAlgebra
import           Grenade.Utils.ListStore

import           Lens.Micro

data BatchNormTape :: Nat -- The number of channels of the tensor
                   -> Nat -- The number of rows of the tensor
                   -> Nat -- The number of columns of the tensor
                   -> Type where
  -- | there is no back propogation when using network for testing, so
  --   so the tape for automatic differentiation is empty
  TestBatchNormTape  :: ( KnownNat channels
                     , KnownNat rows
                     , KnownNat columns)
                     => BatchNormTape channels rows columns -- ^ empty tape

  -- | prevents having to calculate information twice during forward and backwards pass
  TrainBatchNormTape :: ( KnownNat channels
                     , KnownNat rows
                     , KnownNat columns)
                     => [R (channels * rows * columns)]  -- ^ xnorm
                     -> R channels                       -- ^ std
                     -> R channels                       -- ^ running mean
                     -> R channels                       -- ^ running variance
                     -> BatchNormTape channels rows columns

data BatchNorm :: Nat -- The number of channels of the tensor
               -> Nat -- The number of rows of the tensor
               -> Nat -- The number of columns of the tensor
               -> Nat -- momentum
               -> Type where
  BatchNorm :: ( KnownNat channels
               , KnownNat rows
               , KnownNat columns
               , KnownNat momentum)
            => Bool                                 -- ^ True: training, False: testing
            -> BatchNormParams channels             -- ^ gamma and beta
            -> R channels                           -- ^ running mean
            -> R channels                           -- ^ running variance
            -> RealNum                              -- ^ epsilon (used for numerical stability)
            -> ListStore (BatchNormParams channels) -- ^ momentum store
            -> BatchNorm channels rows columns momentum

data BatchNormParams :: Nat -> Type where
  BatchNormParams :: KnownNat channels
                  => R channels -- gamma
                  -> R channels -- beta
                  -> BatchNormParams channels

data BatchNormGrad :: Nat -- The number of channels of the tensor
                   -> Nat -- The number of rows of the tensor
                   -> Nat -- The number of columns of the tensor
                   -> Type where
  BatchNormGrad :: ( KnownNat channels
                   , KnownNat rows
                   , KnownNat columns)
            => R channels -- ^ running mean
            -> R channels -- ^ running variance
            -> R channels -- ^ derivative of loss wrt gamma
            -> R channels -- ^ derivative of loss wrt beta
            -> BatchNormGrad channels rows columns

-- | NFData instances
instance NFData (BatchNormTape channels rows columns) where
  rnf TestBatchNormTape = ()
  rnf (TrainBatchNormTape xnorm std mean var) = rnf xnorm `seq` rnf std `seq` rnf mean `seq` rnf var

instance NFData (BatchNormParams flattenSize) where
  rnf (BatchNormParams gamma beta)
    = rnf gamma `seq` rnf beta

instance NFData (BatchNorm channels rows columns momentum) where
  rnf (BatchNorm training bnparams mean var eps store)
    = rnf training `seq` rnf bnparams `seq` rnf mean `seq` rnf var `seq` rnf eps `seq` rnf store

-- | Show instance

instance Show (BatchNorm channels rows columns momentum) where
  show _ = "Batch Normalization"

-- Serialize instances
instance (KnownNat channels, KnownNat rows, KnownNat columns, KnownNat momentum)
  => Serialize (BatchNorm channels rows columns momentum) where

  put (BatchNorm training bnparams mean var ε store) = do
    put training
    put bnparams
    putListOf put . LA.toList . extract $ mean
    putListOf put . LA.toList . extract $ var
    put ε
    put store

  get = do
    training <- get
    bnparams <- get
    mean     <- maybe (fail "Vector of incorrect size") return . create . LA.fromList =<< getListOf get
    var      <- maybe (fail "Vector of incorrect size") return . create . LA.fromList =<< getListOf get
    ε        <- get
    store    <- get

    return $ BatchNorm training bnparams mean var ε store

instance KnownNat flattenSize => Serialize (BatchNormParams flattenSize) where
  put (BatchNormParams gamma beta) = do
    putListOf put . LA.toList . extract $ gamma
    putListOf put . LA.toList . extract $ beta

  get = do
    gamma <- maybe (fail "Vector of incorrect size") return . create . LA.fromList =<< getListOf get
    beta  <- maybe (fail "Vector of incorrect size") return . create . LA.fromList =<< getListOf get

    return $ BatchNormParams gamma beta

-- | Neural network operations
instance (KnownNat channels, KnownNat rows, KnownNat columns, KnownNat mom, KnownNat (channels * rows * columns))
  => UpdateLayer (BatchNorm channels rows columns mom) where

  type Gradient (BatchNorm channels rows columns mom) = BatchNormGrad channels rows columns
  type MomentumStore (BatchNorm channels rows columns mom) = ListStore (BatchNormParams channels)

  reduceGradient = undefined

  runUpdate opt@OptSGD{} x@(BatchNorm training (BatchNormParams oldGamma oldBeta) _ _ ε store) (BatchNormGrad runningMean runningVar dGamma dBeta)
    = let BatchNormParams oldGammaMomentum oldBetaMomentum = getData opt x store
          VectorResultSGD newGamma newGammaMomentum        = descendVector opt (VectorValuesSGD oldGamma dGamma oldGammaMomentum)
          VectorResultSGD newBeta  newBetaMomentum         = descendVector opt (VectorValuesSGD oldBeta dBeta oldBetaMomentum)
          newStore                                         = setData opt x store (BatchNormParams newGammaMomentum newBetaMomentum)
      in  BatchNorm training (BatchNormParams newGamma newBeta) runningMean runningVar ε newStore

  runUpdate opt@OptAdam{} x@(BatchNorm training (BatchNormParams oldGamma oldBeta) _ _ ε store) (BatchNormGrad runningMean runningVar dGamma dBeta)
    = let [BatchNormParams oldMGamma oldMBeta, BatchNormParams oldVGamma oldVBeta] = getData opt x store
          VectorResultAdam newGamma newMGamma newVGamma                            = descendVector opt (VectorValuesAdam (getStep store) oldGamma dGamma oldMGamma oldVGamma)
          VectorResultAdam newBeta  newMBeta  newVBeta                             = descendVector opt (VectorValuesAdam (getStep store) oldBeta  dBeta  oldMBeta  oldVBeta)
          newStore                                                                 = setData opt x store [BatchNormParams newMGamma newMBeta, BatchNormParams newVGamma newVBeta]
      in  BatchNorm training (BatchNormParams newGamma newBeta) runningMean runningVar ε newStore

  runSettingsUpdate NetworkSettings{trainingActive=training} (BatchNorm _ bnparams mean var ε store) = BatchNorm training bnparams mean var ε store

instance (KnownNat channels, KnownNat rows, KnownNat columns, KnownNat mom)
  => LayerOptimizerData (BatchNorm channels rows columns mom) (Optimizer 'SGD) where

  type MomentumDataType (BatchNorm channels rows columns mom) (Optimizer 'SGD) = BatchNormParams channels
  getData opt x store = head $ getListStore opt x store
  setData opt x store = setListStore opt x store . return
  newData _ _ = BatchNormParams (konst 0) (konst 0)

instance (KnownNat channels, KnownNat rows, KnownNat columns, KnownNat mom)
  => LayerOptimizerData (BatchNorm channels rows columns mom) (Optimizer 'Adam) where

  type MomentumDataType (BatchNorm channels rows columns mom) (Optimizer 'Adam) = BatchNormParams channels
  type MomentumExpOptResult (BatchNorm channels rows columns mom) (Optimizer 'Adam) = [BatchNormParams channels]
  getData     = getListStore
  setData     = setListStore
  newData _ _ = BatchNormParams (konst 0) (konst 0)

instance (KnownNat channels, KnownNat rows, KnownNat columns, KnownNat momentum)
  => RandomLayer (BatchNorm channels rows columns momentum) where
  createRandomWith _ _ = pure initBatchNorm

-- | Initialize a batch norm layer with gamma set to zero for each channel and beta
--   set to one for each channel. Running mean is 0 and running var is 1. The default
--   value for epsilon is 0.00001
initBatchNorm :: forall channels rows columns momentum.
  (KnownNat channels, KnownNat rows, KnownNat columns, KnownNat momentum)
  => BatchNorm channels rows columns momentum
initBatchNorm =
  let ch     = fromIntegral $ natVal (Proxy :: Proxy channels)
      zeroes = replicate ch 0
      ones   = replicate ch 1
      gamma  = vector ones   :: R channels
      beta   = vector zeroes :: R channels
      mean   = vector zeroes :: R channels
      var    = vector ones   :: R channels
      ε      = 0.00001
  in BatchNorm True (BatchNormParams gamma beta) mean var ε mkListStore

instance (KnownNat rows, KnownNat momentum)
  => Layer (BatchNorm 1 1 rows momentum) ('D1 rows) ('D1 rows) where

  type Tape (BatchNorm 1 1 rows momentum) ('D1 rows) ('D1 rows) = BatchNormTape 1 1 rows

  runForwards (BatchNorm True _ _ _ _ _) _
    = error "Cannot train use batch size of 1 with BatchNorm layer during training"

  runForwards (BatchNorm False (BatchNormParams gamma beta) runningMean runningVar ε _) (S1D x)
    = let gamma'       = extract gamma LA.! 0
          beta'        = extract beta LA.! 0
          runningMean' = extract runningMean LA.! 0
          runningVar'  = extract runningVar LA.! 0

          std = sqrt $ runningVar' + ε

          y = dvmap (\x -> ((x - runningMean') / std ) * gamma' + beta') x

      in (TestBatchNormTape, S1D y)

  runBatchForwards bn@(BatchNorm False _ _ _ _ _) xs
    = let outs = map (snd . runForwards bn) xs
      in ([TestBatchNormTape], outs)

  runBatchForwards (BatchNorm True (BatchNormParams gamma beta) runningMean runningVar ε _) xs
    = let [m]             = vectorToList runningMean :: [RealNum]
          [v]             = vectorToList runningVar  :: [RealNum]
          [g]             = vectorToList gamma       :: [RealNum]
          [b]             = vectorToList beta        :: [RealNum]
          mom             = (/ 100) $ fromIntegral $ natVal (Proxy :: Proxy momentum)

          xs'             = map extractV xs

          sample_mean     = batchNormMean xs'                 :: RealNum
          sample_var      = batchNormVariance xs'             :: RealNum

          m'              = mom * m + (1 - mom) * sample_mean :: RealNum
          v'              = mom * v + (1 - mom) * sample_var  :: RealNum
          std             = sqrt $ sample_var + ε             :: RealNum

          x_extracted     = map (\(S1D x) -> x) xs
          x_normalised    = map (dvmap (\a -> (a - sample_mean) / std)) x_extracted
          scaledShifted   = map (dvmap (\a -> g * a + b)) x_normalised
          out             = map S1D scaledShifted

          stdV            = listToVector [std] :: R 1
          runningMeanV    = listToVector [m']  :: R 1
          runningVarV     = listToVector [v']  :: R 1

      in ([TrainBatchNormTape x_normalised stdV runningMeanV runningVarV], out)

  runBatchBackwards (BatchNorm True (BatchNormParams _ _) _ _ _ _) _ _
    = undefined

  runBatchBackwards _ _ _
    = undefined

instance (KnownNat rows, KnownNat columns, KnownNat momentum)
  => Layer (BatchNorm 1 rows columns momentum) ('D2 rows columns) ('D2 rows columns) where

  type Tape (BatchNorm 1 rows columns momentum) ('D2 rows columns) ('D2 rows columns) = BatchNormTape 1 rows columns

  runForwards (BatchNorm True _ _ _ _ _) _
    = error "Cannot train use batch size of 1 with BatchNorm layer during training"

  runForwards (BatchNorm False (BatchNormParams gamma beta) runningMean runningVar ε _) (S2D x)
    = let rows     = fromIntegral $ natVal (Proxy :: Proxy rows)
          columns  = fromIntegral $ natVal (Proxy :: Proxy columns)

          gamma'       = extract gamma
          beta'        = extract beta
          runningMean' = extract runningMean
          runningVar'  = extract runningVar
          mat          = extract x

          y  = batchnorm 1 rows columns ε mat gamma' beta' runningMean' runningVar'
          y' = fromJust . create $ y

      in (TestBatchNormTape, S2D y')

  runBatchForwards bn@(BatchNorm False _ _ _ _ _) xs
    = let outs = map (snd . runForwards bn) xs
      in ([TestBatchNormTape], outs)

  runBatchForwards (BatchNorm True (BatchNormParams gamma beta) runningMean runningVar ε _) xs
    = let [m]             = vectorToList runningMean :: [RealNum]
          [v]             = vectorToList runningVar  :: [RealNum]
          [g]             = vectorToList gamma       :: [RealNum]
          [b]             = vectorToList beta        :: [RealNum]
          mom             = (/ 100) $ fromIntegral $ natVal (Proxy :: Proxy momentum)

          xs'             = map (sflatten . extractM2D) xs

          sample_mean     = batchNormMean xs'                 :: RealNum
          sample_var      = batchNormVariance xs'             :: RealNum

          m'              = mom * m + (1 - mom) * sample_mean :: RealNum
          v'              = mom * v + (1 - mom) * sample_var  :: RealNum
          std             = sqrt $ sample_var + ε             :: RealNum

          x_extracted     = map (\(S2D x) -> x) xs
          x_normalised    = map (dmmap (\a -> (a - sample_mean) / std)) x_extracted
          scaledShifted   = map (dmmap (\a -> g * a + b)) x_normalised
          out             = map S2D scaledShifted

          x_normalised'   = map sflatten x_normalised
          stdV            = listToVector [std] :: R 1
          runningMeanV    = listToVector [m']  :: R 1
          runningVarV     = listToVector [v']  :: R 1

      in ([TrainBatchNormTape x_normalised' stdV runningMeanV runningVarV], out)

  runBatchBackwards (BatchNorm True (BatchNormParams _ _) _ _ _ _) _ _
    = undefined

  runBatchBackwards _ _ _
    = undefined

instance (KnownNat channels, KnownNat rows, KnownNat columns, KnownNat momentum)
  => Layer (BatchNorm channels rows columns momentum) ('D3 rows columns channels) ('D3 rows columns channels) where

  type Tape (BatchNorm channels rows columns momentum) ('D3 rows columns channels) ('D3 rows columns channels) = BatchNormTape channels rows columns

  runForwards (BatchNorm True _ _ _ _ _) _
    = error "Cannot train use batch size of 1 with BatchNorm layer during training"

  runForwards (BatchNorm False (BatchNormParams gamma beta) runningMean runningVar ε _) (S3D x)
    = let rows     = fromIntegral $ natVal (Proxy :: Proxy rows)
          columns  = fromIntegral $ natVal (Proxy :: Proxy columns)
          channels = fromIntegral $ natVal (Proxy :: Proxy channels)

          gamma'       = extract gamma
          beta'        = extract beta
          runningMean' = extract runningMean
          runningVar'  = extract runningVar
          mat          = extract x

          y  = batchnorm channels rows columns ε mat gamma' beta' runningMean' runningVar'
          y' = fromJust . create $ y

      in (TestBatchNormTape, S3D y')

  runBatchForwards bn@(BatchNorm False _ _ _ _ _) xs
    = let outs = map (snd . runForwards bn) xs
      in ([TestBatchNormTape], outs)

  runBatchForwards (BatchNorm True (BatchNormParams gamma beta) runningMean runningVar ε _) xs
    = let ms     = vectorToList runningMean
          vs     = vectorToList runningVar
          gs     = vectorToList gamma
          bs     = vectorToList beta

          cs     = map splitChannels xs :: [[S ('D2 rows columns)]]
          cs'    = transpose cs

          f c g b m v = let gs' = listToVector [g] :: R 1
                            bs' = listToVector [b] :: R 1
                            ms' = listToVector [m] :: R 1
                            vs' = listToVector [v] :: R 1
                            bn' = BatchNorm True (BatchNormParams gs' bs') ms' vs' ε undefined :: BatchNorm 1 rows columns momentum
                        in  runBatchForwards bn' c

          (tapes, outs) = unzip $ zipWith5 f cs' gs bs ms vs
          outs' = transpose outs
      in (combineTapes tapes, map combineChannels outs')
    where
      combineTapes :: [[BatchNormTape 1 i j]] -> [BatchNormTape k i j]
      combineTapes = undefined

  runBatchBackwards (BatchNorm True (BatchNormParams _ _) _ _ _ _) _ _
    = undefined

  runBatchBackwards _ _ _ = undefined

instance OnnxOperator (BatchNorm channels rows columns momentum) where
  onnxOpTypeNames _ = ["BatchNormalization"]

instance (KnownNat channels, KnownNat rows, KnownNat columns, KnownNat momentum) => OnnxLoadable (BatchNorm channels rows columns momentum) where
  loadOnnxNode inits node = case node ^. #input of
    [_, scale, b, mean, var] -> do
      epsilon     <- readFloatAttributeToRealNum "epsilon" node
      loadedScale <- readInitializerVector inits scale
      loadedB     <- readInitializerVector inits b
      loadedMean  <- readInitializerVector inits mean
      loadedVar   <- readInitializerVector inits var

      return $ BatchNorm False (BatchNormParams loadedScale loadedB) loadedMean loadedVar epsilon mkListStore
    _               -> onnxIncorrectNumberOfInputs

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{- | Relevant Vector Machine. Classifier -}
module Numeric.BMML.RVMC
    ( fit
    , predict
    ) where

import           Numeric.LinearAlgebra.Data (Matrix, Vector)
import qualified Numeric.LinearAlgebra.HMatrix as H

data RVMCModel = RVMCModel
    { mean :: Vector Double
    , covariance :: Matrix Double
    }

initAlpha :: Int -> Vector Double
initAlpha = H.konst 1

initY :: Int -> Vector Double
initY = H.konst 0.5

calcB :: Vector Double -> Matrix Double
calcB y =
    H.diag
        (H.cmap
             (\yi ->
                   (sigmoid yi) * (1 - sigmoid yi))
             y)

calcSigma :: Matrix Double -> Matrix Double -> Vector Double -> Matrix Double
calcSigma phi b alpha = H.inv (H.tr phi H.<> b H.<> phi + H.diag alpha)

calcW
    :: Vector Double
    -> Matrix Double
    -> Matrix Double
    -> Matrix Double
    -> Vector Double
    -> Vector Double
    -> Vector Double
calcW mu sigma b phi t y = (sigma H.<> (H.tr phi) H.<> b) H.#> hatt
  where
    hatt = phi H.#> mu + (H.inv b) H.#> (t - y)

calcY :: Vector Double -> Matrix Double -> Vector Double
calcY w phi = H.cmap sigmoid (phi H.#> w)

sigmoid :: Double -> Double
sigmoid x = 1 / (1 + exp (-x))

fit :: Matrix Double -> Vector Double -> RVMCModel
fit phi t =
    let (n,d) = H.size phi
        alpha = initAlpha d
        y = initY n
        b = calcB y
        sigma = calcSigma phi b alpha
    in go
           (2 * floor (sqrt (fromIntegral n)))
           alpha
           y
           b
           (H.konst 0.1 d)
           sigma
  where
    go 0 _ _ _ w sigma =
        RVMCModel
        { mean = w
        , covariance = sigma
        }
    go k alpha y b w sigma = go (k - 1) na ny nb nw ns
      where
        gamma = H.cmap (1 -) (alpha * (H.takeDiag sigma))
        na = gamma / (w * w)
        ny = calcY w phi
        nb = calcB ny
        ns = calcSigma phi nb na
        nw = calcW w sigma b phi t ny

-- | Return probability of phi is from C_1 (t = 1)
predict :: RVMCModel -> Vector Double -> Double
predict (RVMCModel w s) phi = sigmoid ((k sigma_a) * m_a)
  where
    m_a = w `H.dot` phi
    sigma_a = phi `H.dot` (s H.#> phi)
    k x = sqrt (1 + 3 * x / 8)

ext: hs

max_stars_repo_name: DbIHbKA/BMML

max_issues_repo_name: DbIHbKA/BMML

max_forks_repo_name: DbIHbKA/BMML


text:
module Main where

import Graphics.Gloss
import Graphics.Gloss.Data.Color
import Data.Complex
import Data.ByteString (pack)

width = 1200 :: Integer
height = 800 :: Integer

type Pixel = Int
type PixelCoords = (Pixel, Pixel)
type RealCoords = (Double, Double)

type Range a = (a, a)

data State = State {
  xMin :: Float,
  xMax :: Float,
  yMin :: Float,
  yMax :: Float
} deriving (Show)

-- [0, 1024] —> [-1, 1]
mapRange (a, b) (c, d) t = c + ((d - c) / (b - a)) * (t - a)

mapPixelToComplex (xMin, xMax) (yMax, yMin) (n, m) = mapX n :+ mapY m
  where mapX = mapRange (0, fromIntegral width) (xMin, xMax)
        mapY = mapRange (0, fromIntegral height) (yMin, yMax)

window = InWindow "Mandelbrot" (fromIntegral width, fromIntegral height) (0, 0)

backgroundColor = white

fps = 1

initialState = State {
  xMin = -2.25,
  xMax = 1.125,
  yMin = -1.125,
  yMax = 1.125
}

computeNumberOfSteps :: (RealFloat a, Num a, Ord a) => Complex a -> Integer
computeNumberOfSteps z = toInteger $ length $ take 50 iterationsInsideRadius2
  where iterationsInsideRadius2 = takeWhile ((<2) . magnitude) (iterate ((+z) . (**2)) 0)

quadrants (a,b)
  | a == 0 = 1
  | otherwise = 0

-- intsToColors :: [Int] -> [Int]
intsToColors = foldr (\n acc -> (parseColor n) ++ acc) []
  where shade n = ((fromIntegral n) `mod` 10) * 25
        parseColor n = [shade n, shade n, shade n, fromIntegral 255]

steps State { xMin = xMin, xMax = xMax, yMin = yMin, yMax = yMax } = map (computeNumberOfSteps . toComplex) pixels
    where pixels = [(x, y) | y <- [0..(fromIntegral $ height - 1)], x <- [0..(fromIntegral $ width - 1)] ]
          toComplex = mapPixelToComplex (xMin, xMax) (yMax, yMin)

byteString state = pack . intsToColors $ steps state

draw :: State -> Picture
draw state = bitmapOfByteString (fromIntegral width) (fromIntegral height) bitmapFormat (byteString state) True
  where bitmapFormat = BitmapFormat TopToBottom PxRGBA

-- handleInputs = undefined

main :: IO ()
main = display window white $ draw initialState
-- main = play window backgroundColor fps initialState draw (const id) (const id)

max_stars_repo_name: anisometropie/mandelbrot

max_issues_repo_name: anisometropie/mandelbrot

max_forks_repo_name: anisometropie/mandelbrot


text:
{-# LANGUAGE ViewPatterns #-}

module Statistics.Covariance.Shrinkage.MMSE (
    mmseCov
  , shrinkCov
  , sampleCov
  , naiveCov
  , pinCoef
  , lwShrinkageCoef
  , rblwShrinkCovCoef
  , rblwShrinkageCoef
  , oracleShrinkCovCoef
  , oracleShrinkageCoef
  ) where

import Numeric.LinearAlgebra hiding ((<>))
import Numeric.Sum (kbn, sumVector)
import Statistics.Sample (mean)
import qualified Data.Vector.Unboxed as VU

mmseCov :: (Matrix Double -> Double) -> Matrix Double -> Herm Double
mmseCov shrink x = shrinkCov s p f
  where
    s = sampleCov x
    f = naiveCov x
    p = shrink x

shrinkCov :: Herm Double -> Double -> Herm Double -> Herm Double
shrinkCov sample coef target = scale (1 - coef) sample `add` scale coef target

sampleCov :: Matrix Double -> Herm Double
sampleCov = snd . meanCov

avgVarCov :: Herm Double -> Herm Double
avgVarCov (unSym -> s) = trustSym $ scale s' (ident p)
  where
    p = cols s
    s' = mean $ takeDiag s

naiveCov :: Matrix Double -> Herm Double
naiveCov = avgVarCov . sampleCov

pinCoef :: Double -> Double
pinCoef = max 0 . min 1

lwShrinkageCoef :: Matrix Double -> Double
lwShrinkageCoef x = numer / denom
  where
    n = rows x
    p = cols x
    s = unSym $ sampleCov x
    trs2 = trace $ s <> s
    tr2s = sq $ trace s
    frob ix = let c = x ? [ix] in sq . norm_Frob $ unSym (mTm c) - s
    numer = sumVector kbn (VU.generate n frob) / fromIntegral (sq n)
    denom = trs2 - tr2s / fromIntegral p

rblwShrinkCovCoef :: Int -> Herm Double -> Double
rblwShrinkCovCoef (fromIntegral -> n) (unSym -> s) = numer / denom
  where
    p = cols s
    trs2 = trace $ s <> s
    tr2s = sq $ trace s
    numer = (n - 2) / n * trs2 + tr2s
    denom = (n + 2) * (trs2 - tr2s / fromIntegral p)

rblwShrinkageCoef :: Matrix Double -> Double
rblwShrinkageCoef x = rblwShrinkCovCoef n (sampleCov x)
  where n = rows x

data OASIter = OASIter Double (Herm Double)

data OAS = OAS Int Int (Herm Double) (Herm Double)

oracleStep :: OAS -> Herm Double -> OASIter
oracleStep (OAS n p s f) sj = OASIter pj' sj'
  where
    n' = fromIntegral n
    p' = fromIntegral p
    trss = trace $ unSym sj <> unSym s
    tr2s = sq . trace $ unSym sj
    numer = (1 - 2 / p') * trss  + tr2s
    denom = (n' + 1 - 2 / p') * trss + (1 - n' / p') * tr2s
    pj' = numer / denom
    sj' = scale (1 - pj') s `add` scale pj' f

oracleIter :: Double -> OAS -> Herm Double -> Double
oracleIter tol oas = go 1 . oracleStep oas
  where
    go p (OASIter p' sj')
      | abs (p - p') < tol = p'
      | otherwise = go p' $ oracleStep oas sj'

oracleShrinkCovCoef :: Int -> Herm Double -> Double
oracleShrinkCovCoef n s = oracleIter 1e-6 (OAS n p s f) s
  where
    p = cols $ unSym s
    f = avgVarCov s

oracleShrinkageCoef :: Matrix Double -> Double
oracleShrinkageCoef x = oracleShrinkCovCoef n (sampleCov x)
  where n = rows x

sq :: Num a => a -> a
sq a = a * a

trace :: Matrix Double -> Double
trace = sumVector kbn . takeDiag

max_stars_repo_name: tsbattman/coviest

max_issues_repo_name: tsbattman/coviest

max_forks_repo_name: tsbattman/coviest


text:
{-# LANGUAGE DataKinds, FlexibleContexts, TemplateHaskell #-}
module Main where
import qualified Data.Vector.Unboxed as V
import Diagrams.Backend.Rasterific
import Diagrams (dims2D, width, height)
import Frames
import Graphics.Rendering.Chart.Backend.Diagrams (defaultEnv, runBackendR)
import Graphics.Rendering.Chart.Easy
import qualified Pipes as P
import qualified Pipes.Prelude as P
import Statistics.Sample.KernelDensity (kde)

-- Data from http://wwwn.cdc.gov/nchs/nhanes/2005-2006/TRIGLY_D.htm
tableTypes "Trigly" "data/trigly_d.csv"

-- Load the data. Invalid records use zeros as a placeholder.
triglyData :: MonadSafe m => P.Producer Trigly m ()
triglyData = readTable "data/trigly_d.csv" P.>-> P.filter ((> 0) . view lBDLDL)

-- Adapted from a Chart example
fillBetween :: String -> [(a, (b, b))] -> EC l (PlotFillBetween a b)
fillBetween title vs = liftEC $ do
  plot_fillbetween_title .= title
  color <- dissolve 0.5 `fmap` takeColor
  plot_fillbetween_style .= solidFillStyle color
  plot_fillbetween_values .= vs

-- | Plot a semi-transparent KDE with a thick black outline.
mkPlot :: String -> [Int] -> EC (Layout Double Double) ()
mkPlot title xs = do plot (fillBetween title pts)
                     plot . liftEC $ do
                       zoom plot_lines_style $ do
                         line_color .= opaque black
                         line_width *= 2
                       plot_lines_values .= [map (_2 %~ snd) pts]
  where pts = V.toList . uncurry (V.zipWith (\x y -> (x, (0, y)))) . kde 128
            $ V.fromList (map fromIntegral xs)

-- | Plot LDL cholesterol and Triglyceride levels.
mkPlots :: [Record [LBXTR, LBDLDL]] -> EC (Layout Double Double) ()
mkPlots xs = do layout_title .= "Distributions"
                layout_x_axis . laxis_title .= "mg/dL"
                layout_all_font_styles . font_size *= 2
                mkPlot "LDL" ldls
                mkPlot "Triglycerides" tris
  where ldls = map (view lBDLDL) xs
        tris = map (view lBXTR) xs

main :: IO ()
main = do env <- defaultEnv bitmapAlignmentFns 640 480
          let chart2diagram = fst . runBackendR env . toRenderable . execEC
          ldlData <- runSafeT . P.toListM $ triglyData P.>-> P.map rcast
          let d = chart2diagram $ mkPlots ldlData
              sz = dims2D (width d) (height d)
          renderRasterific "plot.png" sz d

max_stars_repo_name: adamConnerSax/Frames

max_issues_repo_name: adamConnerSax/Frames

max_forks_repo_name: adamConnerSax/Frames


text:
module Linalg where

{-# LANGUAGE RankNTypes,
             FlexibleContexts,
             FlexibleInstances,
             MultiParamTypeClasses,
             UndecidableInstances,
             AllowAmbiguousTypes,
             DataKinds,
             ScopedTypeVariables,
#-}

import Data.Matrix
import Data.Bifunctor
import qualified Data.Vector as V
import qualified Data.List.Zipper as ZP
import Data.Complex
import qualified Data.List as DL
import qualified Data.Set as DS
import qualified Data.Map as DM
import GHC.Float
import Data.Maybe
import Data.Euclidean

subtr :: Num a => Matrix a -> Matrix a -> Matrix a
subtr = elementwise (-)

add :: Num a => Matrix a -> Matrix a -> Matrix a
add = elementwise (+)

norm2 :: Floating a => Matrix a -> a
norm2 v = sqrt $ foldr (\x r -> x * x + r) 0 (toList v)

norm2C :: RealFloat a => Matrix (Complex a) -> a
norm2C v = sqrt $ foldr (\x r -> realPart (x * conjugate x) + r) 0 (toList v)

getDiagonal :: Num a => Matrix a -> Matrix a
getDiagonal m = diagonalList (nrows m) 0 $ V.toList $ getDiag m

gershgorinCircles :: Num a => Matrix a -> [(a, a)]
gershgorinCircles m = zip cs rs
    where
        cs   = V.toList $ getDiag m
        diag = getDiagonal m
        rs   = map (foldr (\x r -> abs x + r) 0) (toLists $ m `subtr` diag)

outUnitCircle :: (Num a, Ord a) => Matrix a -> Bool
outUnitCircle m = foldr (\(c, r) acc -> (abs c + r >= 1) || acc) False (gershgorinCircles m)

{- Simple iteration method of solving linear equations systems. Complexity O(n^2) per iteration. -}
simpleIteration :: (Floating a, Ord a) => [[a]] -> [[a]] -> a -> Either Int (Matrix a)
simpleIteration m b = simpleIteration' (fromLists m) (fromLists b)

simpleIteration' :: (Floating a, Ord a) => Matrix a -> Matrix a -> a -> Either Int (Matrix a)
simpleIteration' m b eps = doIterations m' b (zero size 1) eps (outUnitCircle m) 0
    where
          size = nrows m
          m'   = identity size `subtr` m

doIterations m b x eps outCircle cnt
    | outCircle
        && cnt == 20                 = Left 0
    | outCircle
        && norm2 x' >= norm2 x + 1   = doIterations m b x' eps outCircle (cnt + 1)
    | norm2 (x `subtr` x') < eps     = Right x'
    | otherwise                      = doIterations m b x' eps outCircle 0
    where x'                         = (m `multStd` x) `add` b

lowerTriangle :: Num a => Matrix a -> Matrix a
lowerTriangle m = mapPos (\(i, j) e -> if (size - i + 1) + j <= size + 1 then e else 0) m
    where
        size = nrows m

upperTriangle :: Num a => Matrix a -> Matrix a
upperTriangle m = transpose $ m' `subtr` diag
    where
        m'   = lowerTriangle $ transpose m
        diag = getDiagonal m

gaussPartial :: (Floating a, Ord a) => Matrix a -> Matrix a -> Matrix a
gaussPartial m b = fromLists $ map (: []) (foldl gaussStep [] mlist)
    where
        m'                = b <|> m
        mlist             = toLists m'
        gaussStep ans row = ans ++ [xi]
            where
                coef = foldr (\e acc -> if e /= 0 && acc == 0 then e else acc) 0 row
                bi   = head row
                row' = tail row
                sm   = sum $ zipWith (*) ans row'
                xi   = (bi - sm) / coef

{- Gauss-Zeidel method of solving linear equations systems. Complexity O(n^2) per iteration. -}
gaussZeidel :: (Floating a, Ord a) => [[a]] -> [[a]] -> a -> Either Int (Matrix a)
gaussZeidel m b = gaussZeidel' (fromLists m) (fromLists b)

gaussZeidel' :: (Floating a, Ord a) => Matrix a -> Matrix a -> a -> Either Int (Matrix a)
gaussZeidel' m b eps = doGaussZeidel l negu b (zero size 1) eps (outUnitCircle m) 0
    where
        size = nrows m
        l    = lowerTriangle m
        negu = zero size size `subtr` upperTriangle m

doGaussZeidel l negu b x eps outCircle cnt
    | outCircle
        && cnt == 20                            = Left 0
    | outCircle
        && norm2 x' >= norm2 x + 1              = doGaussZeidel l negu b x' eps outCircle (cnt + 1)
    | norm2 ((l `multStd` x) `subtr` b') < eps  = Right x
    | otherwise                                 = doGaussZeidel l negu b x' eps outCircle 0
    where
        b' = (negu `multStd` x) `add` b
        x' = gaussPartial l b'

{- Givens rotation of matrix. Complexity O(n). -}
givensRotation :: Num a => [[a]] -> Int -> Int -> a -> a -> [[a]]
givensRotation m i j c s = zipWith subst [1..] m
    where
        ui    = m !! (i - 1)
        uj    = m !! (j - 1)
        uinew = zipWith (\xi xj -> c * xi + s * xj) ui uj
        ujnew = zipWith (\xi xj -> (- s) * xi + c * xj) ui uj
        subst pos row
            | pos == i  = uinew
            | pos == j  = ujnew
            | otherwise = row

trans :: [[a]] -> [[a]]
trans = toLists . transpose . fromLists

idMatrix :: Num a => Int -> [[a]]
idMatrix sz = toLists $ identity sz

firstNonzero :: (Num a, Ord a) => [a] -> Int -> Int
firstNonzero v k = foldl (\acc (j, x) -> if j >= k && x /= 0 && acc == 0
                                       then j
                                       else acc) 0 (zip [1..] v)

appToPair f (x, y) = (f x, f y)

{- Givens rotation method of QR decomposition. Complexity O(n^3). -}
qrDecompGivens :: (Floating a, Ord a) => [[a]] -> (Matrix a, Matrix a)
qrDecompGivens m = appToPair fromLists $ first trans qr
    where
        idm                = idMatrix $ length m
        qr                 = foldl handler (idm, m) [1..(length m - 1)]
        handler (q, r) k
            | i == 0    = (q, r)
            | otherwise = (givensRotation q'' k i 0 1, givensRotation r'' k i 0 1)
            where
                col        = trans r !! (k - 1)
                i          = firstNonzero col k
                (q'', r'') = fst $ foldl handler' ((q, r), col !! (i - 1)) (zip [1..] col)
                handler' ((q', r'), xi) (j, xj)
                    | j <= i    = ((q', r'), xi)
                    | otherwise = ((givensRotation q' j i c s, givensRotation r' j i c s), (- s) * xj + c * xi)
                        where
                            n = sqrt $ xi * xi + xj * xj
                            c = xi / n
                            s = (- xj) / n

{- Householder matrix multiplication. Complexity O(n^2). -}
multHouseholder :: Num a => [[a]] -> [[a]] -> [[a]]
multHouseholder m' v' = toLists $ m `subtr` prod
    where
        m    = fromLists m'
        v    = fromLists v'
        vt   = transpose v
        prod = scaleMatrix 2 v `multStd` (vt `multStd` m)

{- Householder reflection method of QR decomposition. Complexity O(n^3). -}
qrDecompHouseholder :: (Floating a, Ord a) => [[a]] -> (Matrix a, Matrix a)
qrDecompHouseholder m = appToPair fromLists $ first trans qr
    where
        idm                = idMatrix $ length m
        qr                 = foldl handler (idm, m) [1..(length m - 1)]
        handler (q, r) k
            | norm2 v' == 0 || u == e1 = (q, r)
            | otherwise                = (multHouseholder q v, multHouseholder r v)
            where
                col = trans r !! (k - 1)
                v'  = fromLists $ zipWith (\j x -> if j < k then [0] else [x]) [1..length col] col
                u   = scaleMatrix (1 / norm2 v') v'
                e1  = mapPos (\(j, _) _-> if j == k then 1 else 0) v'
                v   = toLists $ scaleMatrix (1 / norm2 (u `subtr` e1)) (u `subtr` e1)

{- Simple iteration method of calculation the maximum modulo eigenvalue. Complexity O(n^2) per iteration. -}
simpleIterationMaxEV :: RealFloat a => [[Complex a]] -> [[Complex a]] -> a -> Int
                        -> Either Int (Complex a, Matrix (Complex a))
simpleIterationMaxEV m v = doItersMaxEV (fromLists m) (fromLists v)

doItersMaxEV m x eps cnt
    | cnt == 0          = Left 0
    | norm2C diff < eps = Right (ev, x)
    | otherwise         = doItersMaxEV m x' eps (cnt - 1)
    where
        norm = norm2C (m `multStd` x) :+ 0
        x'   = scaleMatrix (1 / norm) (m `multStd` x)
        ev   = head $ head $ toLists $ transpose x `multStd` (m `multStd` x)
        diff = (m `multStd` x) `subtr` scaleMatrix ev x

{- QR-algorithm of calculation the matrix spectrum. Complexity O(n^3) per iteration. -}
qrEV :: (Floating a, Ord a) => [[a]] -> a -> ([a], Matrix a)
qrEV m eps = doItersQrEV m eps (identity $ length m)

doItersQrEV m eps qk
    | lessEps   = (evs, qk)
    | otherwise = doItersQrEV m' eps (qk `multStd` q)
    where
        (q, r)  = qrDecompGivens m
        m'      = toLists $ r `multStd` q
        circles = gershgorinCircles (fromLists m)
        lessEps = foldr (\(_, rd) acc -> (rd < eps) && acc) True circles
        evs     = V.toList $ getDiag $ fromLists m

multHouseholderRight :: Num a => [[a]] -> [[a]] -> [[a]]
multHouseholderRight m' v' = toLists $ m `subtr` prod
    where
        m    = fromLists m'
        v    = fromLists v'
        vt   = transpose v
        prod = scaleMatrix 2 ((m `multStd` v) `multStd` vt)

{- Matrix transformation to tridiagonal form. Complexity O(n^3). -}
getTridiagonal :: (Floating a, Ord a) => [[a]] -> (Matrix a, Matrix a)
getTridiagonal m = appToPair fromLists $ second trans tridiag
    where
        idm                = idMatrix $ length m
        tridiag            = foldl handler (m, idm) [1..(length m - 1)]
        handler (a, q) k
            | norm2 v' == 0 || u == e1 = (a, q)
            | otherwise                = (newa, newq)
            where
                col  = trans a !! (k - 1)
                v'   = fromLists $ zipWith (\j x -> if j <= k then [0] else [x]) [1..length col] col
                u    = scaleMatrix (1 / norm2 v') v'
                e1   = mapPos (\(j, _) _-> if j == k + 1 then 1 else 0) v'
                v    = toLists $ scaleMatrix (1 / norm2 (u `subtr` e1)) (u `subtr` e1)
                newa = multHouseholderRight (multHouseholder a v) v
                newq = multHouseholder q v

cursorp :: ZP.Zipper a -> Int
cursorp (ZP.Zip l _) = length l

rightn :: Int -> ZP.Zipper a -> ZP.Zipper a
rightn 0 z = z
rightn n z = rightn (n - 1) (ZP.right z)

givensRotationZ :: Num a => [ZP.Zipper a] -> Int -> Int -> a -> a -> [ZP.Zipper a]
givensRotationZ m i j c s = zipWith subst [1..] m
    where
        ui      = m !! (i - 1)
        uj      = m !! (j - 1)
        curspUi = cursorp ui
        curspUj = cursorp uj
        uinew   = rightn curspUi $ ZP.fromList $
                   zipWith (\xi xj -> c * xi + s * xj) (ZP.toList ui) (ZP.toList uj)
        ujnew   = rightn curspUj $ ZP.fromList $
                   zipWith (\xi xj -> (- s) * xi + c * xj) (ZP.toList ui) (ZP.toList uj)
        subst pos row
            | pos == i  = uinew
            | pos == j  = ujnew
            | otherwise = row

{- QR decomposition for the tridiagonal matrices. Complexity O(n^2). -}
qrDecompTridiagonal :: (Floating a, Ord a) => [[a]] -> ([(Int, Int, a, a)], Matrix a)
qrDecompTridiagonal m = second (fromLists . fmap ZP.toList) $
                        foldl handler ([], zipped_m) [1..(length m - 1)]
    where
        zipped_m = fmap ZP.fromList m
        handler (gs, r) k
            | n == 0    = (gs, r)
            | otherwise = ((k + 1, k, c, s) : gs, map ZP.right $ givensRotationZ r (k + 1) k c s)
            where
                col = map ZP.cursor r
                xi  = col !! (k - 1)
                xj  = col !! k
                n   = sqrt $ xi * xi + xj * xj
                c   = xi / n
                s   = (- xj) / n

multGivens :: Num a => [(Int, Int, a, a)] -> [[a]] -> [[a]]
multGivens gs m = foldr (\(i, j, c, s) acc -> givensRotation acc i j c s) m gs

{-
QR-algorithm of calculation the matrix spectrum for the tridiagonal matrices.
Complexity O(n^2) per iteration.
-}
qrEVTridiagonal :: (Floating a, Ord a) => [[a]] -> a -> ([a], Matrix a)
qrEVTridiagonal m eps = second (transpose . fromLists) $ doItersQrEVsTridiagonal m eps (idMatrix $ length m)

doItersQrEVsTridiagonal m eps qk
    | lessEps   = (evs, qk)
    | otherwise = doItersQrEVsTridiagonal m' eps q
    where
        (gs, r) = qrDecompTridiagonal m
        m'      = trans $ multGivens gs (trans $ toLists r)
        q       = multGivens gs qk
        circles = gershgorinCircles (fromLists m)
        lessEps = foldr (\(_, rd) acc -> (rd < eps) && acc) True circles
        evs     = V.toList $ getDiag $ fromLists m

doItersQrEVsTridiagonalNTimes m cnt qk
    | cnt == 0  = (evs, qk)
    | otherwise = doItersQrEVsTridiagonalNTimes m' (cnt - 1) q
    where
        (gs, r) = qrDecompTridiagonal m
        m'      = trans $ multGivens gs (trans $ toLists r)
        q       = multGivens gs qk
        evs     = V.toList $ getDiag $ fromLists m

doItersQrMinEVTridiagonal m eps qk
    | radius < eps = (qk, m)
    | otherwise    = doItersQrMinEVTridiagonal m' eps q
    where
        (gs, r) = qrDecompTridiagonal m
        m'      = trans $ multGivens gs (trans $ toLists r)
        q       = multGivens gs qk
        lastRow = last m
        radius  = foldr (\x acc -> abs x + acc) 0 lastRow - abs (last lastRow)

swapMinor :: Num a => Matrix a -> Matrix a -> Matrix a
swapMinor minor m = m' `add` minor'
    where
        minorSz = nrows minor
        mSz     = nrows m
        m'     = mapPos (\(i, j) x -> if i <= minorSz && j <= minorSz
                                  then 0
                                  else x) m
        minor' = extendTo 0 mSz mSz minor

{-
Wilkinson shifts QR-algorithm of calculation the matrix spectrum for the tridiagonal matrices.
Complexity O(n^2) per iteration.
-}
qrEVShifts :: (Floating a, Ord a) => [[a]] -> a -> ([a], Matrix a)
qrEVShifts mt eps = (evs, transpose q)
    where
        sz     = length mt
        (q, d) = foldl handler (identity (length mt), fromLists mt) [0..sz - 2]
        evs    = V.toList $ getDiag d
        handler (qk, m) k
               = (qt `multStd` qk, swapMinor r' m)
               where
                    m'       = submatrix 1 (sz - k) 1 (sz - k) m
                    sz'      = nrows m'
                    lowerSq  = submatrix (sz' - 1) sz' (sz' - 1) sz' m'
                    maxIters = 20
                    s        = head $ tail $ fst $
                               doItersQrEVsTridiagonalNTimes (toLists lowerSq) maxIters (idMatrix 2)
                    m''      = m' `subtr` scaleMatrix s (identity sz')
                    (qt, r)  = first fromLists $ doItersQrMinEVTridiagonal (toLists m'') eps (idMatrix sz)
                    r'       = fromLists r `add` scaleMatrix s (identity sz')

{- Test graphs on nonisomorphism. Complexity O(Time(qrEVShifts)). -}
isomorphic :: (Floating a, Ord a) => [[a]] -> [[a]] -> Bool
isomorphic g1 g2
    | length g1 /= length g2 = False
    | otherwise              = norm2 (g1Spec `subtr` g2Spec) <= 10 * eps
    where
        eps    = 0.00001
        g1'    = toLists $ fst $ getTridiagonal g1
        g2'    = toLists $ fst $ getTridiagonal g2
        g1Spec = fromLists $ map (: []) $ DL.sort $ fst $ qrEVShifts g1' eps
        g2Spec = fromLists $ map (: []) $ DL.sort $ fst $ qrEVShifts g2' eps

getAdjMap :: Num a => Int -> DM.Map Int (DM.Map Int a)
getAdjMap n = adjMap
    where
        buildRow = DM.fromList $ map (\k -> (k, 0)) [0..n - 1]
        adjMap   = DM.fromList $ map (\k -> (k, buildRow)) [0..n - 1]

addEdge (x, y) = DM.update (Just . DM.update (Just . (+ 1)) y) x

buildAdjMatrix :: Num a => Int -> (Int -> DM.Map Int (DM.Map Int a) -> DM.Map Int (DM.Map Int a)) -> Matrix a
buildAdjMatrix n genEdges = m
    where
        adjMap      = getAdjMap n
        adjMap'     = foldr genEdges adjMap [0..n - 1]
        m           = matrix n n fill
        fill (i, j) = x
            where
                i' = i - 1
                j' = j - 1
                x  = fromJust $ DM.lookup (j - 1) $ fromJust $ DM.lookup (i - 1) adjMap'

buildGraph1 :: Int -> Matrix Double
buildGraph1 n = buildAdjMatrix (n * n) genEdges
    where
        genEdges x = foldr (\y acc -> genEdges' (x, y) . acc) id [0..n - 1]
        genEdges' (x, y)
            | x >= n    = id
            | otherwise = addEdge (v, u1) . addEdge (v, u2) .
                          addEdge (v, u3) . addEdge (v, u4) .
                          addEdge (v, u5) . addEdge (v, u6) .
                          addEdge (v, u7) . addEdge (v, u8)
            where
                x1 = (x + 2 * y) `mod` n
                x2 = (x - 2 * y + 3 * n) `mod` n
                x3 = (x + 2 * y + 1) `mod` n
                x4 = (x - 2 * y - 1 + 3 * n) `mod` n
                y1 = (y + 2 * x) `mod` n
                y2 = (y - 2 * x + 3 * n) `mod` n
                y3 = (y + 2 * x + 1) `mod` n
                y4 = (y - 2 * x - 1 + 3 * n) `mod` n
                u1 = x1 * n + y
                u2 = x2 * n + y
                u3 = x3 * n + y
                u4 = x4 * n + y
                u5 = x * n + y1
                u6 = x * n + y2
                u7 = x * n + y3
                u8 = x * n + y4
                v  = x * n + y

inv :: Int -> Int -> Int
inv x md = ((xinv `mod` md) + md) `mod` md
    where
        xinv = snd $ gcdExt x md

buildGraph2 :: Int -> Matrix Double
buildGraph2 p = buildAdjMatrix (p + 1) genEdges
    where
        genEdges x
            | x == p    = addEdge (x, xInv) . addEdge (x, x) . addEdge (x, x)
            | otherwise = addEdge (x, xInv) . addEdge (x, x1) . addEdge (x, x2)
            where
                xInv
                    | x == 0    = p
                    | x == p    = 0
                    | otherwise = inv x p
                x1 = (x + 1) `mod` p
                x2 = (x - 1 + p) `mod` p

{- Calculation the optimal alpha for expander. Complexity O(Time(qrEVShifts)). -}
expanderAlpha :: Int -> Matrix Double -> Double -> Double
expanderAlpha n g d = max (abs ev1) (abs ev2) / d
    where
        eps       = 0.00001
        evs       = fst $ qrEVShifts (toLists $ fst $ getTridiagonal $ toLists g) eps
        sortedEVs = reverse $ DL.sort evs
        ev1       = head $ tail sortedEVs
        ev2       = last sortedEVs

expanderAlpha1 :: Int -> Double
expanderAlpha1 n = expanderAlpha n g 8
    where
        g = buildGraph1 n

expanderAlpha2 :: Int -> Double
expanderAlpha2 n = expanderAlpha (n + 1) g 3
    where
        g = buildGraph2 n


text:
-- the multidimensional minimization example in the GSL manual
import Numeric.GSL
import Numeric.LinearAlgebra
import Graphics.Plot
import Text.Printf(printf)

-- the function to be minimized
f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30

-- exact gradient
df [x,y] = [20*(x-1), 40*(y-2)]

-- a minimization algorithm which does not require the gradient
minimizeS f xi = minimize NMSimplex2 1E-2 100 (replicate (length xi) 1) f xi

-- Numerical estimation of the gradient
gradient f v = [partialDerivative k f v | k <- [0 .. length v -1]]

partialDerivative n f v = fst (derivCentral 0.01 g (v!!n)) where
    g x = f (concat [a,x:b])
    (a,_:b) = splitAt n v

disp = putStrLn . format "  " (printf "%.3f")

allMethods :: (Enum a, Bounded a) => [a]
allMethods = [minBound .. maxBound]

test method = do
    print method
    let (s,p) = minimize method 1E-2 30 [1,1] f [5,7]
    print s
    disp p

testD method = do
    print method
    let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f df [5,7]
    print s
    disp p

testD' method = do
    putStrLn $ show method ++ " with estimated gradient"
    let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f (gradient f) [5,7]
    print s
    disp p

main = do
    mapM_ test [NMSimplex, NMSimplex2]
    mapM_ testD allMethods
    testD' ConjugateFR
    mplot $ drop 3 . toColumns . snd $ minimizeS f [5,7]

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-# LANGUAGE CPP
           , TypeOperators
           , KindSignatures
           , TypeApplications
           , ScopedTypeVariables
           , AllowAmbiguousTypes
           , TypeFamilies
           , FlexibleContexts
           , DeriveFunctor
           , DeriveFoldable
           , DeriveTraversable
           , TemplateHaskell
           , RankNTypes
           , UndecidableInstances
           -- , DeriveAnyClass
           -- , GeneralizedNewtypeDeriving
           -- , DatatypeContexts
  #-}
{-# OPTIONS_GHC -Wall #-}

module Main where

import Control.Applicative
import Control.Monad (unless)
import Data.Complex
import Data.Traversable
import System.Exit (exitFailure)
import Test.QuickCheck (choose, vectorOf)  -- , elements, collect)
import Test.QuickCheck.Arbitrary
import Test.QuickCheck.All (quickCheckAll)

dft :: [Complex Double] -> [Complex Double]
dft xs = [ sum [ x * exp((0.0 :+ (-1.0)) * 2 * pi / lenXs * fromIntegral(k * n))
          | (x, n) <- zip xs [0..]  ]
          | k <- [0..(length xs - 1)]  ]
 where lenXs = fromIntegral $ length xs

type Unop a = a -> a

dft' :: forall f a. ( Applicative f, Traversable f, Sized f
                    , RealFloat a) => Unop (f (Complex a))
dft' xs = omegas (size @f) `cross` xs

cross :: (Functor n, Foldable m, Applicative m, Num a) =>
  n (m a) -> m a -> n a  -- matrix * vector
mat `cross` vec = dot vec <$> mat

dot :: (Foldable f, Applicative f, Num a) =>
  f a -> f a -> a        -- dot product
u `dot` v = sum (liftA2 (*) u v)

newtype (g :. f) a = O (g (f a))
  deriving (Eq, Functor, Foldable, Traversable)

instance (Show a, Show (f a), Show (g String), Functor g) => Show ((g :. f) a) where
  show (O x) = show $ fmap show x

instance (Applicative f, Applicative g) => Applicative (g :. f) where
  pure = O . pure . pure
  O h <*> (O x) = O $ fmap (<*>) h <*> x

class Sized (f :: * -> *) where
    size :: Integer

instance (Sized f, Sized g) => Sized (g :. f) where
    size = size @f * size @g

class FFT f where
  type Reverse f :: * -> *
  fft :: RealFloat a => f (Complex a) -> Reverse f (Complex a)

-- Constraints filled in, by responding to compiler errors incrementally.
instance (Traversable g, Traversable (Reverse g),
          Applicative (Reverse g),
          FFT g,
          Traversable f,
          Applicative f, Applicative (Reverse f),
          FFT f, Sized f, Sized (Reverse g) ) =>
  FFT (g :. f) where
    type Reverse (g :. f) = Reverse f :. Reverse g
    fft = O . fft' . transpose . twiddle . fft' . unO
fft' :: ( Traversable g
        , Applicative (Reverse g)
        , FFT g
        , Traversable f
        , Applicative f
        , RealFloat a ) =>
  g (f (Complex a)) -> Reverse g (f (Complex a))
fft' = transpose . fmap fft . transpose

twiddle :: forall g f a. ( Applicative g, Traversable g, Sized g
                         , Applicative f, Traversable f, Sized f
                         , RealFloat a ) =>
  g (f (Complex a)) -> g (f (Complex a))
twiddle = (liftA2 . liftA2) (*) (omegas (size @(g :. f)))

omegas :: ( Applicative g, Traversable g
          , Applicative f, Traversable f
          , RealFloat a) =>
  Integer -> g (f (Complex a))
-- 'i' was unrecognized.
-- omegas n = powers <$> powers (exp (-i * 2 * pi / fromIntegral n))
omegas n = powers <$> powers (cis (2 * pi / fromIntegral n))

powers :: ( Applicative f, Traversable f
          , Fractional a) => a -> f a
-- I don't have time to digest Conal's "Generic Parallel Scan" talk, right now.
-- powers = fst . lscanAla Product . pure
-- But, I do need to accomodate all applicative functors, as opposed to just lists.
-- So, I'll try the Traversable approach.
powers w = fmap (/ w) . snd . mapAccumL (\x y -> (x * y, x * y)) 1 $ pure w

unO :: (g :. f) a -> g (f a)
unO (O x) = x

-- Copied from Conal's generic-fft library.
transpose :: (Traversable g, Applicative f) => g (f a) -> f (g a)
transpose = sequenceA

-- Quad
newtype Quad a = Quad ((a,a),(a,a))
  deriving (Show, Functor, Foldable, Traversable)
instance Applicative Quad where
  pure x = Quad ((x,x),(x,x))
  Quad ((g,h),(u,v)) <*> Quad ((w,x),(y,z)) = Quad ((g w, h x),(u y, v z))
instance Sized Quad where
  size = 4

-- Pair
data  Pair a = a :# a
  deriving ( Show, Eq, Functor, Foldable, Traversable )
instance Applicative Pair where
  pure x = x :# x
  g :# h <*> (x :# y) = g x :# h y
instance Sized Pair where
  size = 2
instance FFT Pair where
  type Reverse Pair = Pair
  fft (x :# y) = (x + y) :# (x - y)

-- Tree of Pairs
data Tree a = Leaf (Pair a)
            | Branch (Pair (Tree a))
  deriving (Show, Eq, Functor, Foldable, Traversable)

-- Tree building utility.
mkTree :: Integer -> [a] -> Tree a
mkTree n xs | n == 0    = Leaf (head xs :# head (tail xs))
            | otherwise = Branch (mkTree (n-1) (evens xs) :# mkTree (n-1) (odds xs))

-- Pair of Pair building utility.
mkPofP :: [a] -> Pair (Pair a)
-- mkPofP []        = error "Empty list provided!"
-- mkPofP [w]       = error "Wrong number of elements given!"
-- mkPofP [w,x]     = error "Wrong number of elements given!"
-- mkPofP [w,x,y]   = error "Wrong number of elements given!"
-- mkPofP [w,x,y,z] = O ((w :# x) :# (y :# z))
-- mkPofP xs        = O (unO (mkPofP (evens xs)) :# (unO (mkPofP (odds xs))))
mkPofP [w,x,y,z] = (w :# y) :# (x :# z)
-- mkPofP xs        = mkPofP (evens xs) :# mkPofP (odds xs)
mkPofP _         = error "Wrong number of elements given!"

evens :: [a] -> [a]
evens []     = []
evens (x:xs) = x : odds xs

odds :: [a] -> [a]
odds []     = []
odds (_:xs) = evens xs

newtype Id a = Id {unId :: a}
  deriving ( Show, Eq, Functor, Foldable, Traversable )

instance Applicative Id where
  pure = Id
  Id g <*> (Id x) = Id (g x)

instance Sized Id where
  size = 1

instance FFT Id where
  type Reverse Id = Id
  fft (Id x) = Id x

toList :: Foldable f => f a -> [a]
toList = foldr (:) []

assertion1 :: Bool
assertion1 = dft (map (:+ 0.0) [1.0, 0.0, 0.0, 0.0]) == [1.0 :+ 0.0, 1.0 :+ 0.0, 1.0 :+ 0.0, 1.0 :+ 0.0]

assertion2 :: Bool
assertion2 = dft (map (:+ 0.0) [1.0, 0.0, 0.0, 0.0]) == toList (dft' (fmap (:+ 0.0) (Quad ((1.0, 0.0), (0.0, 0.0)))))

assertion3 :: Bool
assertion3 = fft (fmap (:+ (0 :: Double)) (O ((1 :# 0) :# (0 :# 0)))) == dft' (fmap (:+ (0::Double)) (O ((1 :# 0) :# (0 :# 0))))

exeMain :: IO ()
exeMain = do
  print (powers (cis (2 * pi / 4)) :: Quad (Complex Float))
  print assertion1
  print assertion2
  print assertion3

-- QuickCheck types & propositions
newtype FFTTestVal  f g = FFTTestVal {
  getVal :: (g :. f) (Complex Double)
} -- deriving (Show)

instance (Show (f (Complex Double)), Show (g String), Functor g) =>
  Show (FFTTestVal f g) where
    show = show . getVal

instance (Applicative f, Applicative g) => Arbitrary (FFTTestVal f g) where
    arbitrary = do
        depth  <- choose (2, 5)
        -- rvals <- vectorOf (2^depth) (choose (-1.0, 1.0))
        -- ivals <- vectorOf (2^depth) (choose (-1.0, 1.0))
        -- let values = zipWith (:+) rvals ivals
        -- return $ FFTTestVal (fmap (:+ 0.0) (mkP depth))
        -- return $ FFTTestVal (tree2FuncComp $ mkTree depth values)
        return $ FFTTestVal $ O (pure (pure (depth :+ (0.0::Double))))

-- prop_fail = False  -- Test that QuickCheck is, actually, running.

prop_fft_test testVal = fft val == dft' val
  where types = testVal :: FFTTestVal Pair Pair
        val   = getVal testVal

-- Entry point for unit tests.
return []
runTests = $quickCheckAll

testMain :: IO ()
testMain = do
    -- allPass <- $quickCheckAll -- Run QuickCheck on all prop_ functions
    allPass <- runTests
    unless allPass exitFailure
    putStrLn "Success!"

-- This is a clunky, but portable, way to use the same Main module file
-- for both an application and for unit tests.
-- MAIN_FUNCTION is preprocessor macro set to exeMain or testMain.
-- That way we can use the same file for both an application and for tests.
-- #ifndef MAIN_FUNCTION
-- #define MAIN_FUNCTION exeMain
-- #endif
main :: IO ()
-- main = MAIN_FUNCTION
main = testMain


text:
module IQ where

import Data.Complex

type IQ = Complex Float

max_stars_repo_name: hexagonal-sun/ayeQ

max_issues_repo_name: hexagonal-sun/ayeQ

max_forks_repo_name: hexagonal-sun/ayeQ


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

module Trade.Report.SparkLine where

import Data.Time.Clock (diffUTCTime)

import qualified Data.Vector as Vec

import qualified Statistics.Sample as Sample

import qualified Text.Blaze.Html5 as H5
import Text.Blaze.Html5 ((!))
import qualified Text.Blaze.Html5.Attributes as H5A

import Text.Blaze.Svg11 (mkPath, l, m)
import qualified Text.Blaze.Svg11 as S
import qualified Text.Blaze.Svg11.Attributes as A

import Trade.Type.Step.Algorithm

import Trade.Type.DeltaSignal (DeltaSignal(..))
import Trade.Type.Equity (Equity(..))
import Trade.Type.Signal (Timeseries, Signal(..))
import qualified Trade.Type.Signal as Signal

data Config = Config {
  width :: Int
  , height :: Int
  }

defConfig :: Config
defConfig = Config 120 80

type Range = Double -> Double

toRange :: Config -> Timeseries Equity -> Range
toRange conf sig =
  let (_, Equity mi) = Signal.minimum sig
      (_, Equity ma) = Signal.maximum sig
      h = fromIntegral (height conf)
      (lower, upper) = if ma == mi then (0.8, 1.2) else (mi, ma)
      stepH = h / (upper - lower)
  in \y -> (upper-y) * stepH

type Domain = Double -> Double

toDomain :: Config -> [Timeseries Equity] -> (Double, Double, Domain)
toDomain conf sigs =
  let ts = map (\s -> realToFrac (fst (Signal.last s) `diffUTCTime` fst (Signal.head s))) sigs
      tmax = maximum ts
      vs = Vec.fromList ts
      mean = Sample.mean vs
      stdDev = Sample.stdDev vs
      w = fromIntegral (width conf)
  in (mean, stdDev, \x -> x * w / tmax)

svg :: Config -> S.Svg -> S.Svg
svg conf inner =
  S.docTypeSvg
  ! A.version "1.1"
  ! A.width (S.toValue (width conf+1))
  ! A.height (S.toValue (height conf+1))
  $ inner

spark :: Config -> (Double, Double, Domain) -> Timeseries Equity -> S.Svg
spark conf (mean, stdDev, dom) sig@(Signal xs) =
  let t0 = fst (Signal.head sig)
      ran = toRange conf sig
      as = Vec.map (\(t, e) -> (realToFrac (t `diffUTCTime` t0), unEquity e)) xs
      f acc (t, x) =  acc >> l (dom t) (ran x)
      sty = H5A.style (H5.stringValue "stroke:#4444ff;stroke-width:1px;fill:none;")

      coord = do

        S.rect
          ! H5A.style (H5.stringValue "stroke:0px;fill:#ff0000;opacity:0.05;shape-rendering:crispedges;")
          ! A.x (S.toValue (dom (mean - stdDev)))
          ! A.y (S.toValue (0 :: Double))
          ! A.width (S.toValue (dom (2 * stdDev)))
          ! A.height (S.toValue (height conf))

        S.line
          ! H5A.style (H5.stringValue "stroke:#00aa00;shape-rendering:crispedges;")
          ! A.x1 (S.toValue (dom mean))
          ! A.y1 (S.toValue (0 :: Double))
          ! A.x2 (S.toValue (dom mean))
          ! A.y2 (S.toValue (height conf))

        S.line
          ! H5A.style (H5.stringValue "stroke:#000000;shape-rendering:crispedges;")
          ! A.x1 (S.toValue (0 :: Double))
          ! A.y1 (S.toValue (ran 1))
          ! A.x2 (S.toValue (width conf))
          ! A.y2 (S.toValue (ran 1))

        S.line
          ! H5A.style (H5.stringValue "stroke:#000000;shape-rendering:crispedges;")
          ! A.x1 (S.toValue (0 :: Double))
          ! A.y1 (S.toValue (0 :: Double))
          ! A.x2 (S.toValue (0 :: Double))
          ! A.y2 (S.toValue (height conf))

{-
        S.circle
          ! H5A.style (H5.stringValue "fill:#000000;")
          ! A.cx (S.toValue (dom (32914.28571428572)))
          ! A.cy (S.toValue (ran 1))
          ! A.r (S.toValue (5 :: Double))

        (S.text_ (H5.preEscapedToHtml (show lst)))
          ! H5A.style (H5.stringValue "fill:#000000;")
          ! A.x (S.toValue (dom 0))
          ! A.y (S.toValue (ran 1))
          -}

  in svg conf $ do
    coord
    S.path ! sty ! A.d (mkPath (Vec.foldl' f (m 0 (ran 1)) as))

toSparkLine ::
  (Functor f, StepFunction step) =>
  step -> f [DeltaSignal ohlc] -> f [S.Svg]
toSparkLine step mp =
  let eqty = Equity 1
      us = fmap (map (stepFunction step eqty)) mp

      f ss =
        let dom = toDomain defConfig ss
        in map (spark defConfig dom) ss

  in fmap f us

max_stars_repo_name: fphh/trade

max_issues_repo_name: fphh/trade

max_forks_repo_name: fphh/trade


text:
{-# LANGUAGE DataKinds #-}
module CannyEdgeEndPoint where

import           Control.Monad             as M
import           Data.Array.Repa           as R
import           Data.Binary               (decodeFile)
import           Data.ByteString           as BS (readFile)
import           Data.Complex
import           Data.List                 as L
import           DFT.Plan
import           FokkerPlanck.DomainChange
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           GHC.Word
import           Image.IO
import           Image.IO
import           Image.Transform
import           Utils.Parallel
import           OpenCV                    as CV hiding (Z)
import           Data.Array.Unboxed        as AU
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Types
import           Utils.Array

main = do
  args <- getArgs
  let (numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:scale0FreqsStr:scaleFreqsStr:histFilePath:numIterationStr:writeSourceFlagStr:cutoffRadiusEndPointStr:cutoffRadiusStr:reversalFactorStr:inputImgPath:threshold1Str:threshold2Str:pixelDistStr:patchNormFlagStr:patchNormSizeStr:approximatedEigenValueStr:numThreadStr:_) =
        args
      numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scale0Freq = read scale0FreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      scale0Freqs = [-scale0Freq .. scale0Freq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      numIteration = read numIterationStr :: Int
      writeSourceFlag = read writeSourceFlagStr :: Bool
      cutoffRadiusEndPoint = read cutoffRadiusEndPointStr :: Int
      cutoffRadius = read cutoffRadiusStr :: Int
      reversalFactor = read reversalFactorStr :: Double
      threshold1 = read threshold1Str :: Double
      threshold2 = read threshold2Str :: Double
      pixelDist = read pixelDistStr :: Int
      patchNormFlag = read patchNormFlagStr :: Bool
      patchNormSize = read patchNormSizeStr :: Int
      approximatedEigenValue = read approximatedEigenValueStr :: Double
      numThread = read numThreadStr :: Int
      folderPath = "output/test/CannyEdgeEndPoint"
  createDirectoryIfMissing True folderPath
  img <-
    (exceptError . coerceMat . imdecode ImreadUnchanged) <$>
    BS.readFile inputImgPath  :: IO (Mat ('S '[ 'D, 'D]) 'D ('S GHC.Word.Word8))
  let edge =
        exceptError . canny threshold1 threshold2 (Just 7) CannyNormL2 $ img
      edgeRepa = R.map fromIntegral . toRepa $ edge
      (Z :. _ :. rows :. cols) = extent edgeRepa
      nonzeros = L.filter (> 0) . R.toList $ edgeRepa
      avg = L.sum nonzeros / (fromIntegral . L.length $ nonzeros)
      edgeRepaSparse =
        computeS .
        pad [numPoint, numPoint, 1] 0 .
        R.extend (Z :. (1 :: Int) :. All :. All) .
        increasePixelDistance pixelDist .
        R.slice
          (R.map
             (\x ->
                if x >= avg
                  then 255
                  else 0)
             edgeRepa) $
        (Z :. (0 :: Int) :. All :. All)
  plotImageRepa (folderPath </> "edge.png") . ImageRepa 8 . computeS $ edgeRepa
  plotImageRepa (folderPath </> "edge_sparse.png") . ImageRepa 8 $
    edgeRepaSparse
  edgeRepaSparse <-
    (\(ImageRepa _ img) -> img) <$>
    readImageRepa (folderPath </> "edge_sparse.png") False
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z2T0S0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          tao
          theta0Freqs
          thetaFreqs
          scale0Freqs
          scaleFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scale0Freqs
             , L.length theta0Freqs
             , L.length scaleFreqsp
             , L.length thetaFreqs
             ]
             0)
  arrR2Z2T0S0 <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      (cutoff cutoffRadius radialArr)
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  -- planReversal <- makeR2Z2T0S0Plan emptyPlan arrR2Z2T0S0
  -- (plan, pathNormMethod) <-
  --   makePatchNormFilter
  --     planReversal
  --     numPoint
  --     numPoint
  --     patchNormFlag
  --     patchNormSize
  plan <- makeR2Z2T0S0Plan emptyPlan arrR2Z2T0S0
  pathNormMethod <-
    if patchNormFlag
      then do
        let points =
              createIndex2D .
              L.map fst . L.filter (\(_, v) -> v /= 0) . AU.assocs $
              (AU.listArray ((0, 0), (numPoint - 1, numPoint - 1)) . R.toList $
               edgeRepaSparse :: AU.Array (Int, Int) Double)
            ys =
              pointCluster
                (connectionMatrixP
                   (ParallelParams numThread 1)
                   (2 * pixelDist)
                   points) $
              points
        M.zipWithM_
          (\i ->
             plotImageRepa (folderPath </> (printf "Cluster%03d.png" i)) .
             ImageRepa 8)
          [1 :: Int ..] .
          cluster2Array numPoint numPoint $
          ys
        return . PowerMethodConnection $ ys
      else return PowerMethodGlobal
  let endPointFilePath =
        folderPath </>
        (printf
           "EndPoint_%d_%d_%d_%d_%d_%d_%.2f_%.2f_%f.dat"
           numPoint
           (round thetaFreq :: Int)
           (round scaleFreq :: Int)
           (round maxScale :: Int)
           (round tao :: Int)
           cutoffRadiusEndPoint
           thetaSigma
           scaleSigma
           reversalFactor)
      edgeRepaSparseNormalized = normalizeValueRange (0, 1) edgeRepaSparse
      bias =
        R.traverse
          edgeRepaSparseNormalized
          (const
             (Z :. (L.length theta0Freqs) :. (L.length scale0Freqs) :. numPoint :.
              numPoint)) $ \f (Z :. _ :. _ :. i :. j) ->
          f (Z :. (0 :: Int) :. i :. j) :+ 0
      eigenVec =
        R.traverse
          edgeRepaSparseNormalized
          (const
             (Z :. (L.length thetaFreqs) :. (L.length scaleFreqs) :.
              (L.length theta0Freqs) :.
              (L.length scale0Freqs) :.
              numPoint :.
              numPoint)) $ \f (Z :. k :. l :. _ :. _ :. i :. j) ->
          if k == div (L.length thetaFreqs) 2 &&
             l == div (L.length scaleFreqs) 2 &&
             f (Z :. (0 :: Int) :. i :. j) > 0
            then 1 /
                 (fromIntegral $ (L.length theta0Freqs * L.length scale0Freqs)) :+
                 0
            else 0
  -- plotImageRepaComplex (folderPath </> "bias_0.png") .
  --   ImageRepa 8 .
  --   computeS .
  --   R.extend (Z :. (1 :: Int) :. All :. All) .
  --   R.sumS . R.sumS . rotate4D . rotate4D $
  --   bias
  endPointFlag <- doesFileExist endPointFilePath
  endPointR2Z2 <-
    if endPointFlag
      then readRepaArray endPointFilePath
      else (do putStrLn "Couldn't find the endpoint data. Start computing..."
               arrR2Z2T0S0EndPoint <-
                 computeUnboxedP $
                 computeR2Z2T0S0ArrayRadial
                   (cutoff cutoffRadiusEndPoint radialArr)
                   numPoint
                   numPoint
                   1
                   maxScale
                   thetaFreqs
                   scaleFreqs
                   theta0Freqs
                   scale0Freqs
               completionFieldR2Z2 <-
                 computeS <$>
                 powerMethodR2Z2T0S0Reversal
                   plan
                   folderPath
                   numPoint
                   numPoint
                   numOrientation
                   thetaFreqs
                   theta0Freqs
                   numScale
                   scaleFreqs
                   scale0Freqs
                   maxScale
                   arrR2Z2T0S0EndPoint
                   pathNormMethod
                   numIteration
                   writeSourceFlag
                   (printf
                      "_%d_%d_%d_%d_%d_%d_%.2f_%.2f_%f_EndPoint"
                      numPoint
                      (round thetaFreq :: Int)
                      (round scaleFreq :: Int)
                      (round maxScale :: Int)
                      (round tao :: Int)
                      cutoffRadiusEndPoint
                      thetaSigma
                      scaleSigma
                      reversalFactor)
                   0.5
                   reversalFactor
                   approximatedEigenValue
                   (computeS bias)
                   eigenVec
               -- writeRepaArray endPointFilePath completionFieldR2Z2
               return completionFieldR2Z2)
  let endPointR2Z2Biased = R.zipWith (*) endPointR2Z2 bias
      biasMag =
        R.sumS .
        R.sumS .
        rotate4D .
        rotate4D .
        R.map magnitude .
        r2z2Tor2s1rp numOrientation thetaFreqs numScale scaleFreqs $
        endPointR2Z2Biased
  plotImageRepa (folderPath </> "EndPointBias.png") .
    ImageRepa 8 . computeS . R.extend (Z :. (1 :: Int) :. All :. All) $
    biasMag
  -- printf "%f %f\n" reversalFactor (R.sumAllS biasMag)
  -- let rotatedEndBias =
  --       computeUnboxedS $
  --       R.zipWith
  --         (+)
  --         (rotateBiasR2Z2T0S0 90 thetaFreqs endPointR2Z2Biased)
  --         (rotateBiasR2Z2T0S0 (-90) thetaFreqs endPointR2Z2Biased)
  -- powerMethodR2Z2T0S0Bias
  --   plan
  --   folderPath
  --   numPoint
  --   numPoint
  --   numOrientation
  --   thetaFreqs
  --   theta0Freqs
  --   numScale
  --   scaleFreqs
  --   scale0Freqs
  --   arrR2Z2T0S0
  --   PowerMethodNone
  --   numIteration
  --   writeSourceFlag
  --   (printf
  --      "_%d_%d_%d_%d_%d_%d_%.2f_%.2f_%f"
  --      numPoint
  --      (round thetaFreq :: Int)
  --      (round scaleFreq :: Int)
  --      (round maxScale :: Int)
  --      (round tao :: Int)
  --      cutoffRadiusEndPoint
  --      thetaSigma
  --      scaleSigma
  --      reversalFactor)
  --   0.5
  --   rotatedEndBias
  --   eigenVec

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# Language TypeApplications #-}

module LineSpec where

import Test.Hspec
import Test.QuickCheck
import Test.Invariant
import Data.Complex
import Data.Fixed (mod')
import Data.Maybe
import Geometry
import Geometry.Testing

spec :: Spec
spec = describe "Line" $ do
    describe "Affinity" $ do
      it "1" $ (asCmp 1 :: Line) == Line (0,1)
      it "2" $ (asCmp 1 :: Ray) == Ray (0,1)
      it "3" $ (asCmp 1 :: Segment) == Segment (0,1)
      it "4" $ property $ \l -> asCmp (cmp l) # at' (refPoint l) == (l :: Line)
      it "5" $ property $ \l -> asCmp (cmp l) # at' (refPoint l) == (l :: Ray)
      it "5" $ property $ \l -> asCmp (cmp l) # at' (refPoint l) == (l :: Segment)
      it "6" $ property $ \x -> cmp (asCmp x :: Line) == x
      it "7" $ property $ \x -> cmp (asCmp x :: Ray) == x
      it "8" $ property $ \x -> cmp (asCmp x :: Segment) == x

    describe "parametrization" $ do
      it "1" $ property $ \(AnyLine l) x -> (l @-> x) ->@ l ~= x

      it "2" $ property $ \(AnyLine l) x ->
        let p = l @-> x
        in l @-> (p ->@ l) ~= p

      it "3" $ line' @Cmp @Cmp 0 1 @-> 0 == 0
      it "4" $ line' @Cmp @Cmp 0 2 @-> 1 == 2
      it "5" $ line' @Cmp @Cmp 0 2 @-> 0.5 == 1

      it "6" $ property $ \(AnyLine l) x ->
        isJust (l @->? x)

      it "7" $ property $ \(AnyRay r) x ->
        isJust (r @->? x) <==> x >= 0

      it "8" $ property $ \(AnySegment s) x ->
        let x' = 2 * cos x
        in isJust (s @->? x') <==> 0 <= x' && x' <= 1

      it "9" $ property $ \(AnyLine l) (AnyPoint p) -> isJust (p ->@? l)

      it "10" $ property $ \(AnyRay r) (AnyPoint p) ->
        let an = deg . angleValue $ angleBetween r (r # through' p)
        in isJust (p ->@? r)
           <==>
           (0 <= an && an <= 90) || (270 <= an && an <= 360)

      it "11" $ property $ \(AnySegment s) (AnyPoint p) ->
        let (p1, p2) = refPoints s
            an1 = deg . angleValue . innerAngle $ angleWithin p p1 p2
            an2 = deg . angleValue . innerAngle $ angleWithin p p2 p1
        in isJust (p ->@? s)
           <==>
           (0 <= an1 && an1 <= 90) && (0 <= an2 && an2 <= 90)

    describe "containing" $ do
      it "1.1" $
        property $ \(AnyLine l) x -> l `isContaining` (l @-> x)
      it "1.2" $
        property $ \(AnyLine l) -> l `isContaining` refPoint l

      it "2.1" $ aRay `isContaining` ((0,0) :: XY)
      it "2.2" $ aRay `isContaining` ((1,0) :: XY)
      it "2.5" $ not $ aRay `isContaining` ((-1,0) :: XY)
      it "2.6" $ not $ aRay `isContaining` ((1,1) :: XY)

      it "3.1" $ aSegment `isContaining` ((0,0) :: XY)
      it "3.2" $ aSegment `isContaining` ((1,0) :: XY)
      it "3.3" $ aSegment `isContaining` (aSegment @-> 0)
      it "3.4" $ aSegment `isContaining` (aSegment @-> 1)
      it "3.5" $ not $ aSegment `isContaining` ((-1,0) :: XY)
      it "3.6" $ not $ aSegment `isContaining` ((3,0) :: XY)

    describe "tangent and family" $ do
      it "1" $ angle (line' @Cmp @Cmp 0 (1:+1)) == 45
      it "2" $ normal (line' @Cmp @Cmp 0 (1:+2)) 0 == asCmp ((-2):+1)

    describe "unit" $ do
      it "1" $ unit (line' @Cmp @Cmp 0 (3:+4)) == 5
      it "2" $ unit (asCmp 1 :: Ray) == 1
      it "3" $ unit (asCmp 1 :: Segment) == 1
      it "4" $ property $ \l -> unit l == norm (l :: Line)
      it "5" $ property $ \l -> unit l == norm (l :: Ray)
      it "6" $ property $ \l -> unit l == norm (l :: Segment)

    describe "at" $
      it "1" $ aLine # at ((1,2) :: XY) == line (1,2) (2,2)

    describe "along" $ do
      it "1" $ aLine # at (2,3) # along 0 == line (2,3) (3,3)
      it "2" $ aLine # at (2,3) # along 90 == line (2,3) (2,4)
      it "3" $
        property $ \a p l ->
                     let _ = (a :: Direction, p :: Point, l :: Line)
                     in l # at' p # along' a == l # along' a # at' p

    describe "extendToLength" $
      it "1" $ property $ \(AnySegment s) d ->
        unit (s # extendToLength d) ~= abs d

    describe "extendTo" $ do
      it "1" $ property $ \(AnySegment s) (AnyLine l) ->
        (asRay s `isIntersecting` l) ==>
        l `isContaining` ((s # extendTo l) @-> 1)
      it "2" $ property $ \(AnySegment s) (AnyCircle c) ->
        (asRay s `isIntersecting` c) ==>
        c `isContaining` ((s # extendTo c) @-> 1)

    describe "heightTo" $
      it "1" $ property $ \(AnyPoint p) (AnyLine l) ->
        not (l `isContaining` cmp p) ==>
        fromMaybe False $ do
          h <- l # heightFrom p
          return $ l `isContaining` end h && l `isOrthogonal` h

      -- it "2" $ property $ \(AnyPoint p) (AnyCircle c) ->
      --   not (c `isContaining` p) ==>
      --   fromMaybe False $ do
      --     h <- p # heightTo c
      --     let n = normal c (project c (end h))
      --     return $ c `isContaining` (end h) && h `isCollinear` n

    describe "distanceTo" $ do
      it "1" $ origin `distanceTo` aLine == 0
      it "2" $ origin `distanceTo` (aLine # at (0,1)) == 1
      it "3" $ origin `distanceTo` (aLine # at (0,1) # along 45) ~= 1/sqrt 2
      it "4" $ origin `distanceTo` aRay == 0
      it "5" $ origin `distanceTo` (aRay # at (3,4)) == 5
      it "6" $ origin `distanceTo` (aRay # at (-1,0)) == 0
      it "7" $ origin `distanceTo` aSegment == 0
      it "8" $ origin `distanceTo` (aSegment # at (3,4)) ~= 5
      it "9" $ origin `distanceTo` (aSegment # at (-2,0)) ~= 1
      it "10" $ origin `distanceTo` segment (-1,0) (0,1) ~= 1/sqrt 2

    describe "intersections" $ do
      it "1" $
        property $ \(AnyLine l1) (AnyLine l2) ->
             isIntersecting l1 l2
             <==>
             not (l1 `isCollinear` l2) || l1 `isContaining` refPoint l2
      let r = Ray (0, 1:+1)
          s = Segment (0, 1:+1)
      it "2.1" $ r `intersections` (r # translate (3,0)) # null
      it "2.2" $ r `intersections` (r # rotate 90 # translate (2,0) ) == [1:+1]
      it "2.3" $ r `intersections` (r # rotate 90 # translate (-1,0) ) # null
      it "3.1" $ r `intersections` (s # translate (3,0)) # null
      it "3.2" $ r `intersections` (s # rotate 90 # translate (2,0) ) == [1:+1]
      it "3.3" $ r `intersections` (s # rotate 90 # translate (2.1,0) ) # null
      it "3.4" $ r `intersections` (s # rotate 90 # translate (1,1) ) == [1:+1]
      it "4.4" $ s `intersections` (s # rotate 90 # translate (1,1) ) == [1:+1]
      it "4.5" $ s `intersections` (s # rotate 90) == [0]
      it "4.6" $ s `intersections` (s # rotate 90 # translate (1.01,1.01) ) # null
      it "4.7" $ s `intersections` (s # rotate 90 # translate (0.5,0.5) ) ~= [0.5:+0.5]

      it "5.1" $
        property $ \(AnyLine l) (AnyRay r) ->
            if l `isCollinear` r
            then null (l `intersections` r) ||
                 (r `intersections` l ~= [start r]
                   && l `intersections` r ~= [start l])
            else l `intersections` r ~= r `intersections` l

      it "5.2" $
        property $ \(AnyLine l) (AnySegment s) ->
             if l `isCollinear` s
             then null (l `intersections` s) ||
                  (s `intersections` l ~= [start s]
                   && l `intersections` s ~= [start l])
             else l `intersections` s ~= s `intersections` l

      it "5.3" $
          property $ \(AnySegment s) (AnyRay r) ->
             if s `isCollinear` r
             then null (s `intersections` r) ||
                  (r `intersections` s ~= [start r]
                   && s `intersections` r ~= [start s])
             else s `intersections` r ~= r `intersections` s

    describe "clipping" $ do
      it "1" $ aLine # rotate 45 # clipBy aSquare == [Segment (0, 1:+1)]
      it "2" $ aLine # rotate 45 # at (0.5,0.5) # clipBy aSquare ==
                [Segment (0, 1:+1)]
      it "3" $ aRay # rotate 45 # at (0.5,0.5) # clipBy aSquare ==
                [Segment (0.5:+0.5, 1:+1)]
      it "4" $ aSegment # rotate 45 # at (0.1,0.1) # clipBy aSquare ==
                [Segment (0.1:+0.1, 0.8071067811865475:+0.8071067811865475)]
      it "5" $ aLine # rotate 45 # at (-0.5,0) # clipBy aSquare ==
                [Segment (0:+0.5, 0.5:+1)]

max_stars_repo_name: MethaHardworker/geometry

max_issues_repo_name: MethaHardworker/geometry

max_forks_repo_name: MethaHardworker/geometry


text:
{-# LANGUAGE FlexibleContexts #-}
module SentimentAnalysis where

import Numeric.LinearAlgebra as LA
import Data.Bifunctor (Bifunctor(bimap))
import Data.List.Split (splitOn)

hingeLossSingle :: (Ord a, Numeric a) => Vector a -> a -> Vector a -> a -> a
hingeLossSingle x y theta theta0 = max 0 (1 - y * (theta LA.<.> x + theta0))

hingeLossFull :: (Numeric a, Ord a, Num (Vector a), Fractional a) => Matrix a -> Vector a -> Vector a -> a -> a
hingeLossFull x y theta theta0 = loss
    where
        theta0s = size y |> repeat theta0
        ones = size y |> repeat 1.0
        agreement = y * (theta LA.<# LA.tr' x + theta0s)
        losses = LA.cmap (max 0) (ones - agreement)
        loss = sumElements losses / fromIntegral (size losses)

perceptronSingleStepUpdate :: (Ord b, Numeric b, Num (Vector b)) => Vector b -> b -> Vector b -> b -> (Vector b, b)
perceptronSingleStepUpdate x y theta theta0 =
    if y * (theta LA.<.> x + theta0) <= 0
        then (theta + (scalar y * x), theta0 + y)
        else (theta, theta0)

perceptronUpdate :: (Ord t, Numeric t, Num (Vector t)) => [Vector t] -> [t] -> Vector t -> t -> [Int] -> (Vector t, t)
perceptronUpdate _ _ theta theta0 [] = (theta, theta0)
perceptronUpdate xs ys theta theta0 (z:zs) =
    let (theta', theta0') = perceptronSingleStepUpdate (xs !! z) (ys !! z) theta theta0
    in perceptronUpdate xs ys theta' theta0' zs

perceptron :: (Ord t, Numeric t, Num (Vector t)) => Matrix t -> Vector t -> Int -> [Int] -> (Vector t, t)
perceptron x y t indices = iterate step (cols x |> repeat 0, 0) !! t
    where step (theta, theta0) = perceptronUpdate (toRows x) (toList y) theta theta0 indices

averagePerceptronUpdate :: (Ord t, Numeric t, Num (Vector t))
    => [Vector t] -> [t] -> Vector t -> t -> [Int] -> [(Vector t, t)]
averagePerceptronUpdate _ _ theta theta0 [] = []
averagePerceptronUpdate xs ys theta theta0 (z:zs) =
    let (theta', theta0') = perceptronSingleStepUpdate (xs !! z) (ys !! z) theta theta0
    in (theta', theta0') : averagePerceptronUpdate xs ys theta' theta0' zs

averagePerceptronList :: (Ord t, Numeric t, Num (Vector t)) => Matrix t -> Vector t -> [Int] -> [[(Vector t, t)]]
averagePerceptronList x y indices = iterate step [(cols x |> repeat 0, 0)]
    where
        step xs = averagePerceptronUpdate (toRows x) (toList y) theta theta0 indices
            where (theta, theta0) = last xs

averagePerceptron :: (Fractional t, Ord t, Numeric t, Num (Vector t))
                        => Matrix t -> Vector t -> Int -> [Int] -> (Vector t, t)
averagePerceptron x y t indices = (bimap ((/ scalar n) . sum) ((/ n) . sum) . unzip . concat) xss
    where
        n = fromIntegral $ t * rows x
        xss = drop 1 $ take (t + 1) $ averagePerceptronList x y indices

pegasosSingleStepUpdate :: (Ord e, Numeric e, Num (Vector e)) => Vector e -> e -> e -> e -> Vector e -> e -> (Vector e, e)
pegasosSingleStepUpdate x y lambda eta theta theta0 =
    if y * (theta LA.<.> x + theta0) <= 1
        then
            (scalar (1 - eta * lambda) * theta + scalar (eta * y) * x, theta0 + eta * y)
        else
            (scalar (1 - eta * lambda) * theta, theta0)

pegasosUpdate :: (Ord t, Numeric t, Num (Vector t), Floating t)
                    => [Vector t] -> [t] -> t -> t -> Vector t -> t -> [Int] -> (Vector t, t, t)
pegasosUpdate _ _ lambda i theta theta0 [] = (theta, theta0, i)
pegasosUpdate xs ys lambda i theta theta0 (z:zs) =
    let (theta', theta0') = pegasosSingleStepUpdate (xs !! z) (ys !! z) lambda (1 / sqrt i) theta theta0
    in pegasosUpdate xs ys lambda (i + 1) theta' theta0' zs

pegasos :: (Ord t, Numeric t, Floating t, Num (Vector t)) => Matrix t -> Vector t -> Int -> t -> [Int] -> (Vector t, t)
pegasos x y t lambda indices = (theta', theta0')
    where
        (theta', theta0', i') = iterate step (cols x |> repeat 0, 0, 1.0) !! t
        step (theta, theta0, i) = pegasosUpdate (toRows x) (toList y) lambda i theta theta0 indices

main :: IO ()
main = do
    dta <- loadMatrix "toy_data.tsv"
    orderDta <- readFile "200.txt"

    let t = 10
    let lambda = 0.2
    let x = dta ¿ [1, 2]
    let y = flatten $ dta ¿ [0]
    let indices = map read $ splitOn "," orderDta

    let (theta, theta0) = perceptron x y t indices
    putStrLn "Perceptron:"
    print theta
    print theta0

    let (theta, theta0) = averagePerceptron x y t indices
    putStrLn "Average Perceptron:"
    print theta
    print theta0

    let (theta, theta0) = pegasos x y t lambda indices
    putStrLn "Pegasos:"
    print theta
    print theta0


text:
module ReadSpec (spec) where

import           Data.Complex
import           Data.Ratio
import           Read
import           Test.Hspec
import           Types

spec :: Spec
spec =
  let r = extractValue . readExpr in
  describe "readExpr" $ do
    it "parseString" $ do
      r "\"foo\"" `shouldBe` String "foo"
      r "\"foo\\\"bar\\\"baz\"" `shouldBe` String "foo\"bar\"baz"

    it "parseBool" $ do
      r "#t" `shouldBe` Bool True
      r "#f" `shouldBe` Bool False

    it "parseAtom" $
      r "foo" `shouldBe` Atom "foo"

    it "parseHex" $
      r "#xFF" `shouldBe` Number 255

    it "parseDigits" $ do
      r "#o101" `shouldBe` Number 65
      r "#d123" `shouldBe` Number 123
      r "#b101" `shouldBe` Number 5

    it "parseDecimal" $
      r "123" `shouldBe` Number 123

    it "parseCharacter" $ do
      r "#\\newline" `shouldBe` Character '\n'
      r "#\\space" `shouldBe` Character ' '
      r "#\\a" `shouldBe` Character 'a'

    it "parseFloat" $ do
      r "1.0" `shouldBe` Float 1.0
      r "0.1" `shouldBe` Float 0.1
      -- readExpr ".0" `shouldStartWith` String "No match:"

    it "parseRatio" $ do
      r "1/2" `shouldBe` Ratio (1 % 2)
      r "2/4" `shouldBe` Ratio (1 % 2)

    it "parseComplex" $ do
      r "1+2i" `shouldBe` Complex (1.0 :+ 2.0)
      r "1.23+2.46i" `shouldBe` Complex (1.23 :+ 2.46)

    it "parseList" $ do
      r "(a test)" `shouldBe` List [Atom "a", Atom "test"]
      r "(a (nested) test)" `shouldBe` List[Atom "a", List[Atom "nested"], Atom "test"]
      -- readExpr "(a test" `shouldStartWith` "No match:"

    it "parseDottedList" $ do
      r "(a b . c)" `shouldBe` DottedList [Atom "a", Atom "b"] (Atom "c")
      r "(a (b . c) d)" `shouldBe` List [Atom "a", DottedList [Atom "b"] (Atom "c"), Atom "d"]

    it "parseQuoted" $ do
      r "'(a b)" `shouldBe` List [Atom "quote", List [Atom "a", Atom "b"]]
      r "(a 'b c)" `shouldBe` List [Atom "a", List [Atom "quote", Atom "b"], Atom "c"]