text:
{-# LANGUAGE AllowAmbiguousTypes                      #-}
{-# LANGUAGE ApplicativeDo                            #-}
{-# LANGUAGE DeriveDataTypeable                       #-}
{-# LANGUAGE DeriveGeneric                            #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE FlexibleInstances                        #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE KindSignatures                           #-}
{-# LANGUAGE MultiParamTypeClasses                    #-}
{-# LANGUAGE PatternSynonyms                          #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE RecordWildCards                          #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeFamilies                             #-}
{-# LANGUAGE TypeInType                               #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE UndecidableInstances                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

module Backprop.Learn.Model.Regression (
    LinReg(..)
  , LogReg, pattern LogReg
  , LRp(..), lrBeta, lrAlpha, runLRp
  , expandInput, expandOutput, reshapeInput, reshapeOutput
  , ARIMA(..), ARIMAp(..), ARIMAs(..)
  , ARIMAUnroll, arimaUnroll
  , ARIMAUnrollFinal, arimaUnrollFinal
  , AR, MA, ARMA
  ) where

import           Backprop.Learn.Initialize
import           Backprop.Learn.Model.Class
import           Backprop.Learn.Model.Combinator
import           Backprop.Learn.Model.Function
import           Backprop.Learn.Model.Parameter
import           Backprop.Learn.Model.State
import           Control.DeepSeq
import           Control.Monad.Primitive
import           Data.Finite
import           Data.Kind
import           Data.List
import           Data.Maybe
import           Data.Proxy
import           Data.Type.Equality
import           Data.Typeable
import           GHC.Generics                          (Generic)
import           GHC.TypeLits.Compare
import           GHC.TypeLits.Extra
import           GHC.TypeNats
import           Lens.Micro
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop
import           Numeric.LinearAlgebra.Static.Vector
import           Numeric.OneLiner
import           Numeric.Opto.Ref
import           Numeric.Opto.Update hiding            ((<.>))
import           Statistics.Distribution
import           Unsafe.Coerce
import qualified Data.Binary                           as Bi
import qualified Data.Type.Tuple                       as T
import qualified Data.Vector.Generic.Sized             as SVG
import qualified Data.Vector.Sized                     as SV
import qualified Data.Vector.Storable.Sized            as SVS
import qualified Numeric.LinearAlgebra                 as HU
import qualified Numeric.LinearAlgebra.Static          as H
import qualified System.Random.MWC                     as MWC

-- | Multivariate linear regression, from an i-vector to an o-vector.
data LinReg (i :: Nat) (o :: Nat) = LinReg
  deriving Typeable

-- | Mutivariate Logistic regression, from an i-vector to an o-vector.
--
-- Essentially a linear regression postcomposed with the logistic
-- function.
type LogReg i o = RMap (R o) (R o) (LinReg i o)

-- | Constructor for a 'LogReg'
pattern LogReg :: LogReg i o
pattern LogReg <- RM _ LinReg
  where
    LogReg = RM logistic LinReg
{-# COMPLETE LogReg #-}

-- | Linear Regression parameter
data LRp i o = LRp
    { _lrAlpha :: !(R o)
    , _lrBeta  :: !(L o i)
    }
  deriving (Generic, Typeable, Show)

instance NFData (LRp i o)
instance (KnownNat i, KnownNat o) => Initialize (LRp i o)
instance (KnownNat i, KnownNat o) => Additive (LRp i o)
instance (KnownNat i, KnownNat o) => Scaling Double (LRp i o)
instance (KnownNat i, KnownNat o) => Metric Double (LRp i o)
instance (KnownNat i, KnownNat o, Ref m (LRp i o) v) => AdditiveInPlace m v (LRp i o)
instance (KnownNat i, KnownNat o, Ref m (LRp i o) v) => ScalingInPlace m v Double (LRp i o)
instance (KnownNat i, KnownNat o) => Bi.Binary (LRp i o)
instance (KnownNat i, KnownNat o) => Backprop (LRp i o)

lrBeta :: Lens (LRp i o) (LRp i' o) (L o i) (L o i')
lrBeta f lrp = (\w -> lrp { _lrBeta = w }) <$> f (_lrBeta lrp)

lrAlpha :: Lens' (LRp i o) (R o)
lrAlpha f lrp = (\b -> lrp { _lrAlpha = b }) <$> f (_lrAlpha lrp)

runLRp
    :: (KnownNat i, KnownNat o, Reifies s W)
    => BVar s (LRp i o)
    -> BVar s (R i)
    -> BVar s (R o)
runLRp lrp x = (lrp ^^. lrBeta) #> x + (lrp ^^. lrAlpha)

instance (KnownNat i, KnownNat o) => Num (LRp i o) where
    (+)         = gPlus
    (-)         = gMinus
    (*)         = gTimes
    negate      = gNegate
    abs         = gAbs
    signum      = gSignum
    fromInteger = gFromInteger

instance (KnownNat i, KnownNat o) => Fractional (LRp i o) where
    (/)          = gDivide
    recip        = gRecip
    fromRational = gFromRational

instance (KnownNat i, KnownNat o) => Floating (LRp i o) where
    pi    = gPi
    sqrt  = gSqrt
    exp   = gExp
    log   = gLog
    sin   = gSin
    cos   = gCos
    tan   = gTan
    asin  = gAsin
    acos  = gAcos
    atan  = gAtan
    sinh  = gSinh
    cosh  = gCosh
    asinh = gAsinh
    acosh = gAcosh
    atanh = gAtanh

instance (KnownNat i, KnownNat o) => Learn (R i) (R o) (LinReg i o) where
    type LParamMaybe (LinReg i o) = 'Just (LRp i o)

    runLearn _ (J_ p) = stateless (runLRp p)

-- | Adjust an 'LRp' to take extra inputs, initialized randomly.
--
-- Initial contributions to each output is randomized.
expandInput
    :: (PrimMonad m, ContGen d, KnownNat i, KnownNat j, KnownNat o)
    => LRp i o
    -> d
    -> MWC.Gen (PrimState m)
    -> m (LRp (i + j) o)
expandInput LRp{..} d g = LRp _lrAlpha . (_lrBeta H.|||) <$> initialize d g

-- | Adjust an 'LRp' to return extra ouputs, initialized randomly
expandOutput
    :: (PrimMonad m, ContGen d, KnownNat i, KnownNat o, KnownNat p)
    => LRp i o
    -> d
    -> MWC.Gen (PrimState m)
    -> m (LRp i (o + p))
expandOutput LRp{..} d g = do
    newAlpha <- initialize d g
    newBeta  <- initialize d g
    pure (LRp (_lrAlpha H.# newAlpha) (_lrBeta H.=== newBeta))

-- | Premute (or remove) inputs
--
-- Removed inputs will simply have their contributions removed from each
-- output.
reshapeInput
    :: KnownNat i
    => SV.Vector i' (Finite i)
    -> LRp i o
    -> LRp i' o
reshapeInput is p = p { _lrBeta = colsL . fmap (β `SV.index`) $ is }
  where
    β = lCols (_lrBeta p)

-- | Premute (or remove) outputs
reshapeOutput
    :: KnownNat o
    => SV.Vector o' (Finite o)
    -> LRp i o
    -> LRp i o'
reshapeOutput is LRp{..} =
    LRp { _lrAlpha = vecR . SVG.convert . fmap (α `SVS.index`) $ is
        , _lrBeta  = rowsL . fmap (β `SV.index`) $ is
        }
  where
    α = rVec _lrAlpha
    β = lRows _lrBeta

-- | Auto-regressive integrated moving average model.
--
-- ARIMA(p,d,q) is an ARIMA model with p autoregressive history terms on
-- the d-th order differenced history and q error (innovation) history
-- terms.
--
-- It is a @'Learn' Double Double@ instance, and is meant to predict the
-- "next step" of a time sequence.
--
-- In this state, it is a runnable stateful model.  To train, use with
-- 'ARIMAUnroll' to unroll and fix the initial state.
data ARIMA :: Nat -> Nat -> Nat -> Type where
    ARIMA :: { _arimaGenYHist :: forall m. PrimMonad m => MWC.Gen (PrimState m) -> m Double
             , _arimaGenEHist :: forall m. PrimMonad m => MWC.Gen (PrimState m) -> m Double
             }
          -> ARIMA p d q

-- | The "unrolled" and "destated" 'ARIMA' model, which takes a vector of
-- sequential inputs and outputs a vector of the model's expected next
-- steps.
--
-- Useful for actually training 'ARIMA' using gradient descent.
--
-- This /fixes/ the initial error history to be zero (or a fixed stochastic
-- sample), and treats the initial output history to be a /learned
-- parameter/.
--
-- @
-- instance 'Learn' ('SV.Vector' n 'Double') (SV.Vector n Double) ('ARIMAUnroll' p q) where
--     -- | Initial state is a parameter, but initial error history is fixed
--     type 'LParamMaybe' (ARIMAUnroll p q) = 'Just (T2 (ARIMAp p q) (ARIMAs p q))
--     type 'LStateMaybe' (ARIMAUnroll p q) = 'Nothing
-- @
type ARIMAUnroll p d q = DeParamAt (T.T2 (ARIMAp p q) (ARIMAs p d q))
                                   (T.T2 (ARIMAp p q) (T.T2 Double (R (p + d))))
                                   (R q)
                                   (UnrollTrainState (Max (p + d) q) (ARIMA p d q))

-- | 'ARIMAUnroll', but only looking at the final output after running the
-- model on all inputs.
--
-- @
-- instance 'Learn' ('SV.Vector' n 'Double') Double ('ARIMAUnrollFinal' p q) where
--     -- | Initial state is a parameter, but initial error history is fixed
--     type 'LParamMaybe' (ARIMAUnrollFinal p q) = 'Just (T2 (ARIMAp p q) (ARIMAs p q))
--     type 'LStateMaybe' (ARIMAUnrollFinal p q) = 'Nothing
-- @
type ARIMAUnrollFinal p d q = DeParamAt (T.T2 (ARIMAp p q) (ARIMAs p d q))
                                        (T.T2 (ARIMAp p q) (T.T2 Double (R (p + d))))
                                        (R q)
                                        (UnrollFinalTrainState (Max (p + d) q) (ARIMA p d q))

splitHist
    :: T.T2 (ARIMAp p q) (ARIMAs p d q)
    -> (T.T2 (ARIMAp p q) (T.T2 Double (R (p + d))), R q)
splitHist (T.T2 p ARIMAs{..}) = (T.T2 p (T.T2 _arimaYPred _arimaYHist), _arimaEHist)

joinHist
    :: T.T2 (ARIMAp p q) (T.T2 Double (R (p + d)))
    -> R q
    -> T.T2 (ARIMAp p q) (ARIMAs p d q)
joinHist (T.T2 p (T.T2 _arimaYPred _arimaYHist)) _arimaEHist = T.T2 p ARIMAs{..}

-- | Constructor for 'ARIMAUnroll'
arimaUnroll
    :: KnownNat q
    => ARIMA p d q
    -> ARIMAUnroll p d q
arimaUnroll a@ARIMA{..} = DPA
    { _dpaSplit      = splitHist
    , _dpaJoin       = joinHist
    , _dpaParam      = 0
    , _dpaParamStoch = fmap vecR . SVS.replicateM . _arimaGenEHist
    , _dpaLearn      = UnrollTrainState a
    }

-- | Constructor for 'ARIMAUnrollFinal'
arimaUnrollFinal
    :: KnownNat q
    => ARIMA p d q
    -> ARIMAUnrollFinal p d q
arimaUnrollFinal a@ARIMA{..} = DPA
    { _dpaSplit      = splitHist
    , _dpaJoin       = joinHist
    , _dpaParam      = 0
    , _dpaParamStoch = fmap vecR . SVS.replicateM . _arimaGenEHist
    , _dpaLearn      = UnrollFinalTrainState a
    }

-- | 'ARIMA' parmaeters
data ARIMAp :: Nat -> Nat -> Type where
    ARIMAp :: { _arimaPhi      :: !(R p)
              , _arimaTheta    :: !(R q)
              , _arimaConstant :: !Double
              }
           -> ARIMAp p q
  deriving (Generic, Show)

-- | 'ARIMA' state
data ARIMAs :: Nat -> Nat -> Nat -> Type where
    ARIMAs :: { _arimaYPred :: !Double
              , _arimaYHist :: !(R (p + d))
              , _arimaEHist :: !(R q)
              }
          -> ARIMAs p d q
  deriving (Generic, Show)

instance (KnownNat p, KnownNat d, KnownNat q) => Learn Double Double (ARIMA p d q) where
    type LParamMaybe (ARIMA p d q) = 'Just (ARIMAp p q)
    type LStateMaybe (ARIMA p d q) = 'Just (ARIMAs p d q)

    runLearn ARIMA{..} (J_ p) x (J_ s) = (y, J_ s')
      where
        d :: L p (p + d)
        d  = difference
        e  = x - (s ^^. arimaYPred)
        y  = (p ^^. arimaConstant)
           + (p ^^. arimaPhi  ) <.> (constVar d #> (s ^^. arimaYHist))
           + (p ^^. arimaTheta) <.> (s ^^. arimaEHist)
        yHist' = case Proxy @1 %<=? Proxy @(p + d) of
          LE Refl -> single y # constVar dropLast #> (s ^^. arimaYHist)
          NLE _ _ -> 0
        eHist' = case Proxy @1 %<=? Proxy @q of
          LE Refl -> single e # constVar dropLast #> (s ^^. arimaEHist)
          NLE _ _ -> 0
        s' = isoVar3 ARIMAs (\(ARIMAs pr yh eh) -> (pr,yh,eh))
                y
                yHist'
                eHist'

monosquare :: forall n. (n <=? (n ^ 2)) :~: 'True
monosquare = unsafeCoerce Refl

dropLast :: forall n. (KnownNat n, 1 <= n) => L (n - 1) n
dropLast = case monosquare @n of
    Refl -> vecL . SVS.generate $ \ij ->
      let i :: Finite n
          j :: Finite (n - 1)
          (i, j) = separateProduct ij
      in  if fromIntegral @_ @Int i == fromIntegral j
            then 1
            else 0

single :: Reifies s W => BVar s Double -> BVar s (R 1)
single = konst

difference'
    :: Int                  -- ^ initial
    -> Int                  -- ^ target
    -> HU.Matrix Double     -- ^ target x initial
difference' n m = foldl' go (HU.ident m) [m + 1 .. n]
  where
    go x k = x HU.<> d k
    d k = HU.build (k-1, k) $ \i j ->
        case round @_ @Int (j - i) of
          0 -> 1
          1 -> -1
          _ -> 0

difference :: forall n m. (KnownNat n, KnownNat m) => L n (n + m)
difference = fromJust . H.create $ difference' (n + m) n
  where
    n = fromIntegral $ natVal (Proxy @n)
    m = fromIntegral $ natVal (Proxy @m)

instance NFData (ARIMAp p q)
instance NFData (ARIMAs p d q)

instance Num (ARIMAp p q) where
    (+)         = gPlus
    (-)         = gMinus
    (*)         = gTimes
    negate      = gNegate
    abs         = gAbs
    signum      = gSignum
    fromInteger = gFromInteger

instance Num (ARIMAs p d q) where
    (+)         = gPlus
    (-)         = gMinus
    (*)         = gTimes
    negate      = gNegate
    abs         = gAbs
    signum      = gSignum
    fromInteger = gFromInteger

instance Fractional (ARIMAp p q) where
    (/)          = gDivide
    recip        = gRecip
    fromRational = gFromRational

instance Fractional (ARIMAs p d q) where
    (/)          = gDivide
    recip        = gRecip
    fromRational = gFromRational

instance Floating (ARIMAp p q) where
    pi    = gPi
    sqrt  = gSqrt
    exp   = gExp
    log   = gLog
    sin   = gSin
    cos   = gCos
    asin  = gAsin
    acos  = gAcos
    atan  = gAtan
    sinh  = gSinh
    cosh  = gCosh
    asinh = gAsinh
    acosh = gAcosh
    atanh = gAtanh

instance Floating (ARIMAs p d q) where
    pi    = gPi
    sqrt  = gSqrt
    exp   = gExp
    log   = gLog
    sin   = gSin
    cos   = gCos
    asin  = gAsin
    acos  = gAcos
    atan  = gAtan
    sinh  = gSinh
    cosh  = gCosh
    asinh = gAsinh
    acosh = gAcosh
    atanh = gAtanh

instance Additive (ARIMAp p q) where
    (.+.)   = gAdd
    addZero = gAddZero
instance Additive (ARIMAs p d q) where
    (.+.)   = gAdd
    addZero = gAddZero

instance (KnownNat p, KnownNat q) => Scaling Double (ARIMAp p q)
instance (KnownNat p, KnownNat d, KnownNat q) => Scaling Double (ARIMAs p d q)

instance (KnownNat p, KnownNat q) => Metric  Double (ARIMAp p q)
instance (KnownNat p, KnownNat d, KnownNat q) => Metric  Double (ARIMAs p d q)

instance (KnownNat p, KnownNat q, Ref m (ARIMAp p q) v) => AdditiveInPlace m v (ARIMAp p q)
instance (KnownNat p, KnownNat d, KnownNat q, Ref m (ARIMAs p d q) v) => AdditiveInPlace m v (ARIMAs p d q)

instance (KnownNat p, KnownNat q, Ref m (ARIMAp p q) v) => ScalingInPlace m v Double (ARIMAp p q)
instance (KnownNat p, KnownNat d, KnownNat q, Ref m (ARIMAs p d q) v) => ScalingInPlace m v Double (ARIMAs p d q)

instance (KnownNat p, KnownNat q) => Initialize (ARIMAp p q)
instance (KnownNat p, KnownNat d, KnownNat q) => Initialize (ARIMAs p d q)

instance (KnownNat p, KnownNat q) => Bi.Binary (ARIMAp p q)
instance (KnownNat p, KnownNat d, KnownNat q) => Bi.Binary (ARIMAs p d q)

instance (KnownNat p, KnownNat q) => Backprop (ARIMAp p q)
instance (KnownNat p, KnownNat d, KnownNat q) => Backprop (ARIMAs p d q)

arimaPhi :: Lens (ARIMAp p q) (ARIMAp p' q) (R p) (R p')
arimaPhi f a = (\x' -> a { _arimaPhi = x' } ) <$> f (_arimaPhi a)

arimaTheta :: Lens (ARIMAp p q) (ARIMAp p q') (R q) (R q')
arimaTheta f a = (\x' -> a { _arimaTheta = x' } ) <$> f (_arimaTheta a)

arimaConstant :: Lens' (ARIMAp p q) Double
arimaConstant f a = (\x' -> a { _arimaConstant = x' } ) <$> f (_arimaConstant a)

arimaYPred :: Lens' (ARIMAs p d q) Double
arimaYPred f a = (\x' -> a { _arimaYPred = x' } ) <$> f (_arimaYPred a)

arimaYHist :: Lens' (ARIMAs p d q) (R (p + d))
arimaYHist f a = (\x' -> a { _arimaYHist = x' } ) <$> f (_arimaYHist a)

arimaEHist :: Lens (ARIMAs p d q) (ARIMAs p d q') (R q) (R q')
arimaEHist f a = (\x' -> a { _arimaEHist = x' } ) <$> f (_arimaEHist a)

-- | Autoregressive model
type AR p = ARIMA p 0 0

-- | Moving average model
type MA = ARIMA 0 0

-- | Autoregressive Moving average model
type ARMA p = ARIMA p 0

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
module Clustering.Measures where

import           Data.Maybe ( fromMaybe )
import qualified Control.Foldl as L

import qualified Data.Vector as V
import Numeric.LinearAlgebra

import Clustering.Hierarchical
import Clustering.Data

meanCD :: Fractional c => c -> [c] -> c
meanCD a = L.fold L.mean . map (abs . (a -))

maxCD :: (Ord c, Fractional c) => c -> [c] -> c
maxCD a = fromMaybe 0 . L.fold L.maximum . map (abs . (a -))

getPatternMatrix :: [V.Vector Float] -> Matrix Float
getPatternMatrix [] = undefined
getPatternMatrix xs@(x:_) = (length xs >< length x) $  concatMap V.toList xs

xs :: [V.Vector Float]
xs = map V.fromList [[2, 5],
                     [6, 4],
                     [5, 3],
                     [2, 2],
                     [1, 4],
                     [5, 2],
                     [3, 3],
                     [2, 3]]

xs2 :: [V.Vector Float]
xs2 = map V.fromList [[1, 1],
                      [2, 1],
                      [5, 4],
                      [6, 5],
                      [6.5, 6]]

euc :: V.Vector Float -> V.Vector Float -> Float
-- euc v1 v2 = sqrt $ V.sum $ V.zipWith (\a b ->  (a - b)^2) v1 v2
euc = ((sqrt . V.sum) .) . V.zipWith (\a b ->  (a - b)^2)

meanCDeuc :: V.Vector Float -> [V.Vector Float] -> Float
meanCDeuc v l = euc v $ ccentre l

ccentre :: (Fractional b, Foldable t) => t (V.Vector b) -> V.Vector b
ccentre l = V.map (/(fromIntegral $ length l))
            $ foldl1 (V.zipWith (+)) l

meanCCD :: [V.Vector Float] -> [V.Vector Float] -> Float
meanCCD l1 l2 = euc (ccentre l1) (ccentre l2)

meanAgg :: Fractional c => (t -> a -> c) -> t -> [a] -> c
meanAgg f a = L.fold L.mean . map (f a)

reprDD :: Dendrogram (V.Vector Float )
    -> Dendrogram (V.Vector Float )  -> Float
reprDD DNil DNil = 0
reprDD DNil _ = 0
reprDD _ DNil = 0
reprDD DNode{repr=r1} DNode{repr=r2} = euc r1 r2

meanDRU :: Dendrogram (V.Vector Float )
    -> Dendrogram (V.Vector Float)  -> V.Vector Float
meanDRU DNil DNil = undefined
meanDRU DNil _ = undefined
meanDRU _ DNil = undefined
meanDRU DNode{repr=r1} DNode{repr=r2} = ccentre [r1, r2]

ext: hs

max_stars_repo_name: zitkat/clustering

max_issues_repo_name: zitkat/clustering

max_forks_repo_name: zitkat/clustering


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Statistics.PCA
-- Copyright   :  (c) A. V. H. McPhail 2010, 2014, 2017
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- Principal Components Analysis
--
-----------------------------------------------------------------------------

module Numeric.Statistics.PCA (
                               pca, pcaN, pcaTransform, pcaReduce, pcaReduceN
                          ) where

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

import qualified Data.Array.IArray as I

import Prelude hiding ((<>))
import Data.List(sortBy)
import Data.Ord(comparing)
import Numeric.LinearAlgebra
import Numeric.GSL.Statistics

--import Numeric.Statistics

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

-- | find the principal components of multidimensional data greater than
--    the threshhold
pca :: I.Array Int (Vector Double)     -- the data
    -> Double                         -- eigenvalue threshold
    -> (Vector Double, Matrix Double) -- Eignevalues, Principal components
pca d q = let d' = fmap (\x -> x - (scalar $ mean x)) d -- remove the mean from each dimension
              d'' = fromColumns $ I.elems d'
              (_,vec',uni') = svd d''
              vec = toList vec'
              uni = toColumns uni'
              v' = zip vec uni
              v = filter (\(x,_) -> x > q) v'  -- keep only eigens > than parameter
              (eigs,vs) = unzip v
          in (fromList eigs,fromColumns vs)

-- | find N greatest principal components of multidimensional data
--    according to size of the eigenvalue
pcaN :: I.Array Int (Vector Double)    -- the data
     -> Int                            -- number of components to return
     -> (Vector Double, Matrix Double) -- Eignevalues, Principal components
pcaN d n = let d' = fmap (\x -> x - (scalar $ mean x)) d -- remove the mean from each dimension
               d'' = fromColumns $ I.elems d'
               (_,vec',uni') = svd d''
               vec = toList vec'
               uni = toColumns uni'
               v' = zip vec uni
               v = take n $ reverse $ sortBy (comparing fst) v'
               (eigs,vs) = unzip v
           in (fromList eigs,fromColumns vs)

v1 = fromList [1,2,3,4,5,6::Double]
v2 = fromList [2,3,4,5,6,7::Double]
v3 = fromList [3,4,5,6,7,8::Double]

a = fromColumns [v1,v2,v3]
b = I.listArray (1,3::Int) [v1,v2,v3] :: I.Array Int (Vector Double)

-- | perform a PCA transform of the original data (remove mean)
-- |     Final = M^T Data^T
pcaTransform :: I.Array Int (Vector Double)    -- ^ the data
             -> Matrix Double                  -- ^ the principal components
             -> I.Array Int (Vector Double)    -- ^ the transformed data
pcaTransform d m = let d' = fmap (\x -> x - (scalar $ mean x)) d -- remove the mean from each dimension
                   in I.listArray (1,cols m) $ toRows $ (tr' m) <> (fromRows $ I.elems d')

-- | perform a dimension-reducing PCA modification,
--     using an eigenvalue threshhold
pcaReduce :: I.Array Int (Vector Double)      -- ^ the data
          -> Double                           -- ^ eigenvalue threshold
          -> I.Array Int (Vector Double)      -- ^ the reduced data
pcaReduce d q = let u = fmap (scalar . mean) d
                    d' = zipWith (-) (I.elems d) (I.elems u)
                    d'' = fromColumns d'
                    (_,vec',uni') = svd d''
                    vec = toList vec'
                    uni = toColumns uni'
                    v' = zip vec uni
                    v = filter (\(x,_) -> x > q) v'  -- keep only eigens > than parameter
                    m = fromColumns $ snd $ unzip v
                 in I.listArray (I.bounds d) $ zipWith (+) (toRows $ m <> (tr' m) <> fromRows d') (I.elems u)

-- | perform a dimension-reducing PCA modification, using N components
pcaReduceN :: I.Array Int (Vector Double)      -- ^ the data
           -> Int                              -- ^ N, the number of components
           -> I.Array Int (Vector Double)      -- ^ the reduced data, with n principal components
pcaReduceN d n = let u = fmap (scalar . mean) d
                     d' = zipWith (-) (I.elems d) (I.elems u)
                     d'' = fromColumns d'
                     (_,vec',uni') = svd d''
                     vec = toList vec'
                     uni = toColumns uni'
                     v' = zip vec uni
                     v = take n $ reverse $ sortBy (comparing fst) v'
                     m = fromColumns $ snd $ unzip v
                  in I.listArray (I.bounds d) $ zipWith (+) (toRows $ m <> (tr' m) <> fromRows d') (I.elems u)

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

max_stars_repo_name: amcphail/hstatistics

max_issues_repo_name: amcphail/hstatistics

max_forks_repo_name: amcphail/hstatistics


text:
-- TODO: Avoid recomputing norms and reg const

{-# LANGUAGE BangPatterns, FlexibleContexts, NamedFieldPuns #-}

module School.Unit.LogSoftMax
( logSoftMax ) where

import Numeric.LinearAlgebra (Container, Matrix, R, Vector, add, atIndex,
                              build, cmap, maxElement, size, sumElements)
import School.Unit.Unit (Unit(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitParams (UnitParams(..))
import School.Unit.UnitGradient (UnitGradient(..))
import School.Utils.LinearAlgebra (mapRows, sumCols)

logSoftMax :: Unit R
logSoftMax = Unit
  { apply = lsmApply
  , deriv = lsmDeriv
  }

getNorm :: (Container Vector a, Floating a)
        => a
        -> Vector a
        -> a
getNorm reg = log
      . sumElements
      . cmap (exp . flip (-) reg)

handleRow :: (Container Vector a, Floating a)
          => Vector a
          -> Vector a
handleRow v =
  let !reg = maxElement v
      !norm = getNorm reg v
  in cmap (flip (-) (norm + reg)) v

lsmApply :: UnitParams R
         -> UnitActivation R
         -> UnitActivation R
lsmApply EmptyParams (BatchActivation inMatrix) =
  BatchActivation $ mapRows handleRow inMatrix
lsmApply _ (ApplyFail msg) = ApplyFail msg
lsmApply _ _ = ApplyFail "Failure when applying log softmax unit"

failure :: String
        -> ( UnitGradient R
           , UnitParams R
           )
failure msg = (GradientFail msg, EmptyParams)

builder :: (Container Vector a, RealFrac a)
        => Matrix a
        -> Vector a
        -> Vector a
        -> (a -> a -> a)
builder exps norms gradCols dj dk =
  let j = round dj :: Int
      k = round dk :: Int
  in (-1) * atIndex exps (j, k)
   * atIndex gradCols j
   / atIndex norms j

lsmDeriv :: UnitParams R
         -> UnitGradient R
         -> UnitActivation R
         -> ( UnitGradient R
            , UnitParams R
            )
lsmDeriv EmptyParams
         (BatchGradient inGrad)
         (BatchActivation input) =
  (outGrad, EmptyParams) where
    !reg = maxElement input
    !exps = cmap (exp . flip (-) reg) input
    !norms = sumCols exps
    !gradCols = sumCols inGrad
    delta = build (size input) (builder exps norms gradCols)
    outGrad = BatchGradient $ add inGrad delta
lsmDeriv _ (GradientFail msg) _ = failure msg
lsmDeriv _ _ (ApplyFail msg) = failure msg
lsmDeriv _ _ _ = failure "Failure when differentiating log softmax unit"

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
{-# LANGUAGE DeriveGeneric              #-}
{-# LANGUAGE FlexibleContexts           #-}
{-# LANGUAGE FlexibleInstances          #-}
{-# LANGUAGE GADTs                      #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE InstanceSigs               #-}
{-# LANGUAGE KindSignatures             #-}
{-# LANGUAGE LambdaCase                 #-}
{-# LANGUAGE MultiParamTypeClasses      #-}
{-# LANGUAGE RankNTypes                 #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE StandaloneDeriving         #-}
{-# LANGUAGE TypeFamilies               #-}
{-# LANGUAGE TypeInType                 #-}
{-# LANGUAGE TypeOperators              #-}
{-# LANGUAGE UndecidableInstances       #-}

module TensorOps.Backend.NTensor
  ( NTensor
  , NTensorL
  , NTensorV
  )
  where

import           Control.DeepSeq
import           Control.Monad.Primitive
import           Data.Distributive
import           Data.Kind
import           Data.List.Util
import           Data.Nested
import           Data.Proxy
import           Data.Singletons
import           Data.Singletons.Prelude.List hiding (Length, Reverse)
import           Data.Type.Combinator
import           Data.Type.Combinator.Util
import           Data.Type.Length                    as TCL
import           Data.Type.Length.Util               as TCL
import           Data.Type.Product                   as TCP
import           Data.Type.Product.Util              as TCP
import           Data.Type.Sing
import           Data.Type.Uniform
import           GHC.Generics
import           Statistics.Distribution
import           System.Random.MWC
import           TensorOps.NatKind
import           TensorOps.Types
import           Type.Class.Higher
import           Type.Class.Higher.Util
import           Type.Class.Witness
import           Type.Family.List
import qualified Data.Type.Vector                    as TCV
import qualified Data.Type.Vector.Util               as TCV
import qualified Data.Vector.Sized                   as VS
import qualified Type.Family.Nat                     as TCN

newtype NTensor :: (k -> Type -> Type) -> Type -> [k] -> Type where
    NTensor
        :: { getNVec :: Nested v js a }
        -> NTensor v a js

deriving instance (Num a, SingI ns, Nesting1 Proxy Functor v, Nesting1 Sing Applicative v)
        => Num (NTensor v a ns)
deriving instance Generic (NTensor v a ns)
deriving instance (NFData a, Nesting Proxy NFData v) => NFData (NTensor v a ns)
instance (NFData a, Nesting Proxy NFData v) => Nesting1 w NFData (NTensor v a) where
    nesting1 _ = Wit

instance (NFData a, Nesting Proxy NFData v) => NFData1 (NTensor v a)

genNTensor
    :: forall k (v :: k -> Type -> Type) (ns :: [k]) (a :: Type). Vec v
    => Sing ns
    -> (Prod (IndexN k) ns -> a)
    -> NTensor v a ns
genNTensor s f = NTensor $ genNested s f
{-# INLINE genNTensor #-}

genNTensorA
    :: forall k (v :: k -> Type -> Type) (ns :: [k]) (a :: Type) f. (Applicative f, Vec v)
    => Sing ns
    -> (Prod (IndexN k) ns -> f a)
    -> f (NTensor v a ns)
genNTensorA s f = NTensor <$> genNestedA s f
{-# INLINE genNTensorA #-}

indexNTensor
    :: forall k (v :: k -> Type -> Type) (ns :: [k]) (a :: Type). Vec v
    => Prod (IndexN k) ns
    -> NTensor v a ns
    -> a
indexNTensor i = indexNested i . getNVec
{-# INLINE indexNTensor #-}

overNVec2
    :: (Nested v ns a -> Nested w ms b -> Nested u os c)
    -> NTensor v a ns
    -> NTensor w b ms
    -> NTensor u c os
overNVec2 f x y = NTensor $ f (getNVec x) (getNVec y)
{-# INLINE overNVec2 #-}

ntNVec
    :: Functor f
    => (Nested v ns a -> f (Nested w ms b))
    -> NTensor v a ns
    -> f (NTensor w b ms)
ntNVec f = fmap NTensor . f . getNVec
{-# INLINE ntNVec #-}

nvecNT
    :: Functor f
    => (NTensor v a ns -> f (NTensor w b ms))
    -> Nested v ns a
    -> f (Nested w ms b)
nvecNT f = fmap getNVec . f . NTensor
{-# INLINE nvecNT #-}

overNVec
    :: (Nested v ns a -> Nested v ms a)
    -> NTensor v a ns
    -> NTensor v a ms
overNVec f = getI . ntNVec (I . f)
{-# INLINE overNVec #-}

underNVec
    :: (NTensor v a ns -> NTensor v a ms)
    -> Nested v ns a
    -> Nested v ms a
underNVec f = getNVec . f . NTensor

instance
      ( Vec (v :: k -> Type -> Type)
      , RealFloat a
      , Nesting1 Proxy Functor      v
      , Nesting1 Sing  Applicative  v
      , Nesting1 Proxy Foldable     v
      , Nesting1 Proxy Traversable  v
      , Nesting1 Sing  Distributive v
      , Eq1 (IndexN k)
      ) => Tensor (NTensor v a) where
    type ElemT  (NTensor v a) = a

    liftT
        :: forall (n :: TCN.N) (o :: [k]). SingI o
        => (TCV.Vec n a -> a)
        -> TCV.Vec n (NTensor v a o)
        -> NTensor v a o
    liftT f = NTensor . liftNested f . fmap getNVec
    {-# INLINE liftT #-}

    sumT = sum'
    {-# INLINE sumT #-}

    scaleT α = overNVec (fmap (α *))
    {-# INLINE scaleT #-}

    transp
        :: forall ns. (SingI ns, SingI (Reverse ns))
        => NTensor v a ns
        -> NTensor v a (Reverse ns)
    -- transp = overNVec (transpose sing)
    transp = overNVec (transpose' (singLength sing) sing)
    {-# INLINE transp #-}

    gmul
        :: forall ms os ns. SingI (Reverse os ++ ns)
        => Length ms
        -> Length os
        -> Length ns
        -> NTensor v a (ms         ++ os)
        -> NTensor v a (Reverse os ++ ns)
        -> NTensor v a (ms         ++ ns)
    gmul lM lO lN = overNVec2 (gmul' lM lO lN)
                      \\ sO'
                      \\ sN
      where
        sO' :: Sing (Reverse os)
        sN  :: Sing ns
        (sO', sN) = splitSing (TCL.reverse' lO)
                              (sing :: Sing (Reverse os ++ ns))
    {-# INLINE gmul #-}

    -- TODO: this is a dumb implementation
    diag
        :: forall n ns. SingI (n ': ns)
        => Uniform n ns
        -> NTensor v a '[n]
        -> NTensor v a (n ': ns)
    diag u d
        = genNTensor sing (\i -> case TCV.uniformVec (prodToVec I (US u) i) of
                                   Nothing     -> 0
                                   Just (I i') -> indexNTensor (i' :< Ø) d
                          )
            \\ (produceEq1 :: Eq1 (IndexN k) :- Eq (IndexN k n))
    {-# INLINE diag #-}

    -- TODO: this can be just done using genNTensor
    getDiag
        :: forall n ns. SingI '[n]
        => Uniform n ns
        -> NTensor v a (n ': n ': ns)
        -> NTensor v a '[n]
    getDiag u = overNVec (diagNV sing u)
                  \\ sHead (sing :: Sing '[n])
    {-# INLINE getDiag #-}

    genRand
        :: forall m d (ns :: [k]). (ContGen d, PrimMonad m, SingI ns)
        => d
        -> Gen (PrimState m)
        -> m (NTensor v a ns)
    genRand d g = generateA (\_ -> realToFrac <$> genContVar d g)
    {-# INLINE genRand #-}

    generateA
        :: forall f ns. (Applicative f, SingI ns)
        => (Prod (IndexN k) ns -> f a)
        -> f (NTensor v a ns)
    generateA = genNTensorA sing
    {-# INLINE generateA #-}

    (!) = flip indexNTensor
    {-# INLINE (!) #-}

    ixRows
        :: Applicative f
        => Length ms
        -> Length os
        -> (Prod (IndexN k) ms -> NTensor v a ns -> f (NTensor v a os))
        -> NTensor v a (ms ++ ns)
        -> f (NTensor v a (ms ++ os))
    ixRows l _ f = ntNVec $ fmap joinNested . nIxRows l (\i -> nvecNT (f i))
    {-# INLINE ixRows #-}

    sumRows
        :: forall n ns. SingI ns
        => NTensor v a (n ': ns)
        -> NTensor v a ns
    sumRows = overNVec sumRowsNested
                \\ (nesting1 Proxy :: Wit (Foldable (v n)))
    {-# INLINE sumRows #-}

    mapRows
        :: Length ns
        -> (NTensor v a ms -> NTensor v a ms)
        -> NTensor v a (ns ++ ms)
        -> NTensor v a (ns ++ ms)
    mapRows l f = overNVec $ joinNested . mapNVecSlices (underNVec f) l
    {-# INLINE mapRows #-}

type NTensorL = NTensor (Flip2 TCV.VecT   I) Double
type NTensorV = NTensor (Flip2 VS.VectorT I) Double

max_stars_repo_name: mstksg/tensor-ops

max_issues_repo_name: mstksg/tensor-ops

max_forks_repo_name: mstksg/tensor-ops


text:
{-# LANGUAGE FlexibleContexts #-} {- Container (Vector a) => ..., etc. -}
{-# LANGUAGE RankNTypes       #-} {- For type traces w/ ST             -}
{-# LANGUAGE Strict           #-} {- Lazy evaluation is no good here   -}
{-# LANGUAGE TemplateHaskell  #-} {- Since we're using lenses aplenty  -}

{-
   TODO LIST:
   * Make sure basic functionality works, write some tests.
   * Use Numeric.LinearAlgebra.Devel to guarantee that ALL matrix operations
     are zero-copy or low allocation overhead.
   * Leverage chunk-based parallelism (lots of embarassing parallelism!)
-}

module Optimization.NelderMead
  ( nelderMeadAt
  , nelderMeadFull
  , nelderMead
  ) where

import Control.Applicative
import Control.Lens
import Control.Monad
import Control.Monad.ST

import Data.Maybe
import Data.Ord   ( comparing )
import Data.STRef ( STRef
                  , modifySTRef
                  , newSTRef
                  , readSTRef
                  , writeSTRef )

import Data.Vector.Storable          as VS  hiding ( take
                                                   , forM_
                                                   , (!) )
import Data.Vector.Algorithms.Intro  ( sortBy )
import Data.Vector.Mutable           as VM  hiding ( take )
import Data.Vector.Storable.Mutable  as VSM hiding ( take )

import Numeric.LinearAlgebra.Data    as LA
import Numeric.LinearAlgebra.HMatrix as HM
import Numeric.LinearAlgebra.Devel   as MutM

import Optimization.Types

-- One point of the simplex and its associated score.
data Vertex s = Vertex {
  which    :: Int,
  score    :: s
  }

-- Ranked vertices for NM search.
data VertexRank a b = VertexRank {
  bestOf         :: Vertex b,
  nextBestOf     :: Vertex b,
  worstOf        :: Vertex b,
  goodCentroidOf :: Vector a
  }

-- Create a VertexRank from the available data.
pickBestAndWorstVerts ::
  (Numeric a, Fractional a, Ord b) =>
  --(Container Vector a, Element a, Ord b, Transposable (Simplex a) (Simplex a)) =>
  (Vector a -> b)  -> -- A function from a point on the simplex to some orderable type.
  Simplex a        -> -- The simplex to consider.
  VertexRank a b      -- Vertex rank.
pickBestAndWorstVerts f s = runST $ do
  let
    nDims  = rows s
    points = tr s

  scores <- VM.new nDims
  forM_ [0..(nDims - 1)] (\i -> VM.unsafeWrite scores i (i, f $ points ! i))
  sortBy (comparing snd) scores
  let
    mvertex j = liftM (uncurry Vertex) (VM.unsafeRead scores j)
  h  <- mvertex 0
  sh <- mvertex 1
  l  <- mvertex (nDims - 1)

  return $ VertexRank h sh l (center $ allButIth (which l) s)

-- Increment a termination criterion by one iteration.
increment :: TerminationCriterion a -> TerminationCriterion a
increment p
  | isJust (view maxit p) = over maxit (fmap pred) p
  | otherwise             = p

-- Compute the center of the NM Simplex.
center ::
  (Numeric a, Fractional a) =>
  Simplex a ->
  Vector a
center s = runST $ do
  centroid <- VSM.new (rows s)

  forM_ [0..((rows s) - 1)] (\i -> do
    forM_ [0..((cols s) - 1)] (\j -> do
      let
        u = fromIntegral (cols s)
      VSM.unsafeWrite centroid i ((s `atIndex` (i,j)) / u) ))

  freeze centroid

-- Create a simplex with all points but the i-th.
allButIth :: (Numeric a) => Int -> Simplex a -> Simplex a
allButIth l s
  | l == 0        = subMatrix (0,1) (nDims, nDims) s
  | l == (cols s) = subMatrix (0,0) (nDims, nDims) s
  | otherwise     = ls ||| rs where
      nDims = rows s
      ls = subMatrix (0, l - 1) (nDims, l - 1)     s
      rs = subMatrix (0, l + 1) (nDims, nDims - l) s

-- The centroid of all points except whichever is worst.
goodCentroid :: (Numeric a, Fractional a, Ord b) =>
  VertexRank a b ->
  Simplex a      ->
  Vector a
goodCentroid r splx = center $ allButIth ((which . worstOf) r) splx

-- Get the i-th point of a simplex.
nthPoint :: (Numeric a, Element a) => Int -> Simplex a -> Vector a
nthPoint i s = flatten $ subMatrix (0, i) (rows s, 1) s

-- Get an estimate of simplex size by computing the sum of all
-- elements in the simplex, translated so that its centroid
-- lies at zero.
simplexSize :: (Numeric a, Fractional a) =>
  Simplex a ->
  a
simplexSize s = fastSimplexL1 (moveToOrigin s) where
  fastSimplexL1 m = runST $ do
    acc <- newSTRef 0
    mut <- unsafeThawMatrix m
    forM_ [0..((rows m) - 1)] $ \i -> do
      forM_ [0..((cols m) - 1)] $ \j -> do
        elem <- unsafeReadMatrix mut i j
        modifySTRef acc (+ elem)
    readSTRef acc

moveToOrigin :: (Numeric a, Fractional a) =>
  Simplex a ->
  Simplex a
moveToOrigin s = runSTMatrix $ do
  sm  <- thawMatrix s
  ctr <- unsafeThawVector (center s)

  forM_ [0..((rows s) - 1)] $ \i -> do
    c_entry <- unsafeReadVector ctr i
    forM_ [0..((cols s) - 1)] $ \j -> do
      m_entry <- unsafeReadMatrix sm i j
      unsafeWriteMatrix sm i j (m_entry - c_entry)

  return sm

-- Check if current simplex state and termination criterion indicate that we should
-- be finished.
isDone tc simplex = noMoreIters || closeEnough || neitherValid where
  noMoreIters  = (fromMaybe 1 (view maxit tc)) == 0
  thresh       = view threshold tc
  closeEnough  = fromMaybe False (fmap ((simplexSize simplex) <) thresh)
  neitherValid = (isNothing (view maxit tc)) && (isNothing (view threshold tc))

-- Update the i-th point of simplex sx with vector p.
updateSimplex :: (Numeric a) => Int -> Vector a -> Simplex a -> Simplex a
updateSimplex i v sx = runSTMatrix $ do
  sxm <- thawMatrix sx
  vm  <- unsafeThawVector v
  forM_ [0..((cols sx) - 1)] $ \d -> do
    updated <- unsafeReadVector vm d
    unsafeWriteMatrix sxm d i updated
  return sxm

-- Reflect the simplex's "worst" point around the hyperplane defined by the
-- other points. If this constitutes a significant improvement, then consider
-- "expanding" even farther in the same direction. Otherwise, give up.
reflectAndExpand :: forall a b s. (Num a, Numeric a, Fractional a, Ord b) =>
  NMStepParameters a              ->
  (Vector a -> b)                 ->
  STRef s (Simplex a)             ->
  STRef s (VertexRank a b)        ->
  ST    s Bool
reflectAndExpand params f s vrank = do
  r  <- readSTRef vrank
  sx <- readSTRef s

  let
    x0      = goodCentroid r sx
    x_best  = nthPoint ((which . bestOf) r) sx
    sx_best = (score . bestOf)     r -- Score of best x
    sx_next = (score . nextBestOf) r -- Score of next best x

    a   = view alpha params
    xr  = x0 + ((x0 - x_best) `scaleBy` a) -- x reflected
    sxr = f xr                             -- score of x reflected

  -- How does the reflected point fare?
  if (sxr < sx_next && sxr > sx_best)
    then do
    modifySTRef s (updateSimplex ((which . worstOf) r) xr)
    return True
    else do

    -- Consider expansion.
    if (sxr < sx_best)
      then do
      let
        g   = view gamma params
        xe  = x0 + ((xr - x0) `scaleBy` g)
        sxe = f xe

      -- Is expansion better?
      if (sxe < sxr)
        then do
        modifySTRef s (updateSimplex ((which . worstOf) r) xe)
        return True
        else do
        modifySTRef s (updateSimplex ((which . worstOf) r) xr)
        return True
      else do
      return False

-- Move the worst point of the simplex inward toward the centroid of all the
-- better points. If this doesn't result in an improvement, don't do anything.
contract :: forall a b s. (Numeric a, Fractional a, Ord b) =>
  NMStepParameters a       ->
  (Vector a -> b)          ->
  STRef s (Simplex a)      ->
  STRef s (VertexRank a b) ->
  ST s Bool
contract params f s vrank = do
  r  <- readSTRef vrank
  sx <- readSTRef s
  let
    x0       = goodCentroid r sx
    x_worst  = nthPoint ((which . worstOf) r) sx
    sx_worst = (score . worstOf) r

    rh  = view rho params
    xc  = x0 + ((x_worst - x0) `scaleBy` rh)
    sxc = f xc

  if (sxc < sx_worst)
    then do
    modifySTRef s (updateSimplex ((which . worstOf) r) xc)
    return True
    else do
    return False

-- Move all points of a simplex slightly towards the best point.
shrink :: forall a b s. (Num a, Fractional a, Storable a, Ord b) =>
  NMStepParameters a       ->
  (Vector a -> b)          ->
  STRef s (Simplex a)      ->
  STRef s (VertexRank a b) ->
  ST s Bool
shrink params f s vrank = do
  sx <- readSTRef s
  r  <- readSTRef vrank
  let
    sig    = view sigma r
    x_best = nthPoint ((which . bestOf) r) sx

  forM_ (\i -> do
            modifySTRef s (\simplex ->
                             if (i == ((which . bestOf) r))
                             then
                               simplex
                             else
                               let
                                 xi  = nthPoint i simplex
                                 xi' = x_best + (sig * (xi - x_best)) in
                                 updateSimplex i xi' simplex
                          )
        ) [0..((cols sx) - 1)]
  return True

-- |Predict the location of a minimum of an objective function using the Nelder-Mead simplex
--  search. This function permits tuning of all parameters, including the starting simplex.
nelderMeadFull :: (Num a, Fractional a, Storable a, Ord b) =>
  NMStepParameters a     -> -- ^Nelder-Mead step parameters
  TerminationCriterion a -> -- ^Reason we would terminate
  (Vector a -> b)        -> -- ^Objective function
  Simplex a              -> -- ^Starting simplex
  Vector a                  -- ^Predicted minimum
nelderMeadFull params term_conditions func initial = (center . runST) $ do

  simplex  <- newSTRef initial
  vtx_rank <- newSTRef (pickBestAndWorstVerts simplex)

  nelderMeadStep params term_conditions func simplex vtx_rank where

    nelderMeadStep p tc f s vrank = do
      -- Update rank and see if we're done.
      modifySTRef vrank (pickBestAndWorstVerts func)
      finished <- isDone <$> readSTRef s

      if (finished)
        then do
        rnk <- readSTRef vrank
        sx  <- readSTRef s
        return $ nthPoint ((which . bestOf) rnk) sx

        -- Attempt to improve simplex.
        else do
        reflected <- reflectAndExpand p f s vrank
        unless reflected $ do
          contracted <- contract p f s vrank
          unless contracted $ do
            shrink p f s vrank

        nelderMeadStep p (iterate tc) f s vrank

-- |Predict the location of a minimum of an objective function using the
--  Nelder-Mead simplex search. This version is like nelderMeadFull, but accepts
--  a starting point rather than an entire simplex.
nelderMeadAt :: (Num a, Fractional a, Storable a, Ord b) =>
  NMStepParameters a     -> -- ^Nelder-Mead step parameters
  TerminationCriterion a -> -- ^Reason we would terminate
  (Vector a -> b)        -> -- ^Objective function
  Vector a               -> -- ^Starting point.
  Vector a                  -- ^Predicted minimum.
nelderMeadAt p t f i = nelderMeadFull p t f (makeSimplexAround i e) where
  e = view epsilon p

-- |Minimize an objective function using the Nelder-Mead simplex search. This
--  invokation uses "default" step parameters.
nelderMead :: (Num a, Fractional a, Storable a, Ord b) =>
  TerminationCriterion a -> -- ^Reason we would terminate
  (Vector a -> b)        -> -- ^Objective function
  Vector a                -- ^Initial guess
nelderMead tc f x = nelderMeadAt defaultStep tc f x

{-
   Construct a simplex around a point by creating another vertex epsilon away from the initial
   guess along each vector of the canonical basis.
-}
makeSimplexAround :: (Element a) =>
  Vector a -> -- One vertex of the simplex
  a        -> -- Epsilon used to expand the simplex along the canonical basis
  Simplex a   -- Resultant simplex.
makeSimplexAround initial = fromBlocks [[asColumn initial, cloned + identity]]
  where
    dimension = size initial
    cloned    = asColumns (take dimension $ repeat initial)
    identity  = ident dimension

scaleBy :: (Num a) => Vector a -> a -> Vector a
scaleBy v a = runST $ do
  vm <- thawVector
  forM_ [0..((VS.length v) - 1)] $ \i -> do
    e <- unsafeReadVector vm i
    unsafeWriteVector vm i (e * a)
  freezeVector vm

max_stars_repo_name: agbrooks/buzzwords

max_issues_repo_name: agbrooks/buzzwords

max_forks_repo_name: agbrooks/buzzwords