text:
{-# LANGUAGE AllowAmbiguousTypes       #-}
{-# LANGUAGE DataKinds                 #-}
{-# LANGUAGE DeriveGeneric             #-}
{-# LANGUAGE PolyKinds                 #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE GADTs                     #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE OverloadedStrings         #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# LANGUAGE TypeApplications          #-}
{-# LANGUAGE TypeFamilies              #-}
{-# LANGUAGE TypeOperators             #-}
{-# LANGUAGE StandaloneDeriving        #-}
{-# LANGUAGE UndecidableInstances        #-}

module BlueRipple.Model.MRP where

import qualified BlueRipple.Data.Keyed         as K
import qualified BlueRipple.Data.DataFrames    as BR
import           BlueRipple.Data.CountFolds
import qualified Control.Foldl                 as FL
import           Control.Monad                  ( join )
import qualified Control.Monad.State           as State
import qualified Data.Array                    as A
import           Data.Function                  ( on )
import qualified Data.List                     as L
import qualified Data.Set                      as S
import qualified Data.IntMap                   as IM
import qualified Data.Map                      as M
import           Data.Maybe                     ( isJust
                                                , fromJust
                                                )
import qualified Data.Text                     as T
import qualified Data.Profunctor               as P
import qualified Data.Serialize                as SE
import           Data.Serialize.Text ()
import qualified Data.Vector                   as V
import qualified Data.Vector.Storable          as VS

import qualified Frames                        as F
import qualified Frames.Melt                   as F
import qualified Frames.InCore                 as FI
import qualified Data.Vinyl                    as V
import qualified Data.Vinyl.TypeLevel          as V

import qualified Control.MapReduce             as MR
import qualified Frames.Transform              as FT
import qualified Frames.Folds                  as FF
import qualified Frames.MapReduce              as FMR
import qualified Frames.Enumerations           as FE
--import qualified Frames.Misc                  as FU
import qualified Frames.Serialize              as FS

import qualified Knit.Report                   as K hiding (elements)
import qualified Polysemy.Error                as P
                                                ( mapError )
import qualified Polysemy                      as P
                                                ( raise )
import           Text.Pandoc.Error             as PE
import qualified Text.Blaze.Colonnade          as BC
import qualified Text.Show

import qualified Numeric.LinearAlgebra         as LA

import qualified Data.IndexedSet               as IS
import qualified Numeric.GLM.ProblemTypes      as GLM
import qualified Numeric.GLM.ModelTypes        as GLM
import qualified Numeric.GLM.FunctionFamily    as GLM
import           Numeric.GLM.MixedModel        as GLM
import qualified Numeric.GLM.Bootstrap         as GLM
import qualified Numeric.GLM.Report            as GLM
import qualified Numeric.GLM.Predict           as GLM
import qualified Numeric.GLM.Confidence        as GLM
import qualified Numeric.SparseDenseConversions
                                               as SD

import qualified Relude.Extra as Relude

import qualified Statistics.Types              as ST
import           GHC.Generics                   ( Generic, Rep )

{- Moved to BlueRipple.Data.CountFolds
-- map reduce folds for counting
type Count = "Count" F.:-> Int
type Successes = "Successes" F.:-> Int
type MeanWeight = "MeanWeight" F.:-> Double
type VarWeight = "VarWeight" F.:-> Double
type WeightedSuccesses = "WeightedSuccesses" F.:-> Double
type UnweightedSuccesses =  "UnweightedSuccesses" F.:-> Int

binomialFold
  :: (F.Record r -> Bool) -> FL.Fold (F.Record r) (F.Record '[Count, Successes])
binomialFold testRow =
  let successesF = FL.premap (\r -> if testRow r then 1 else 0) FL.sum
  in  (\n s -> n F.&: s F.&: V.RNil) <$> FL.length <*> successesF

countFold
  :: forall k r d
   . ( Ord (F.Record k)
     , FI.RecVec (k V.++ '[Count, Successes])
     , k F.⊆ r
     , d F.⊆ r
     )
  => (F.Record d -> Bool)
  -> FL.Fold (F.Record r) [F.FrameRec (k V.++ '[Count, Successes])]
countFold testData = MR.mapReduceFold
  MR.noUnpack
  (FMR.assignKeysAndData @k)
  (FMR.foldAndAddKey $ binomialFold testData)

type CountCols = '[Count, UnweightedSuccesses, WeightedSuccesses, MeanWeight, VarWeight]

zeroCount :: F.Record CountCols
zeroCount = 0 F.&: 0 F.&: 0 F.&: 1 F.&: 0 F.&: V.RNil

weightedBinomialFold
  :: (F.Record r -> Bool)
  -> (F.Record r -> Double)
  -> FL.Fold (F.Record r) (F.Record CountCols)
weightedBinomialFold testRow weightRow =
  let wSuccessesF =
        FL.premap (\r -> if testRow r then weightRow r else 0) FL.sum
      successesF  = FL.premap (\r -> if testRow r then 1 else 0) FL.sum
      meanWeightF = FL.premap weightRow FL.mean
      varWeightF  = FL.premap weightRow FL.variance
      f s ws mw = if mw < 1e-6 then realToFrac s else ws / mw -- if meanweight is 0 but
  in  (\n s ws mw vw -> n F.&: s F.&: f s ws mw F.&: mw F.&: vw F.&: V.RNil)
      <$> FL.length
      <*> successesF
      <*> wSuccessesF
      <*> meanWeightF
      <*> varWeightF

weightedCountFold
  :: forall k r d
   . (Ord (F.Record k)
     , FI.RecVec (k V.++ CountCols)
     , k F.⊆ r
     , k F.⊆ (k V.++ r)
     , d F.⊆ (k V.++ r)
     )
  => (F.Record r -> Bool) -- ^ count this row?
  -> (F.Record d -> Bool) -- ^ success ?
  -> (F.Record d -> Double) -- ^ weight
  -> FL.Fold (F.Record r) (F.FrameRec (k V.++ CountCols))
weightedCountFold {- filterData testData weightData-} = weightedCountFoldGeneral (F.rcast @k)
{-  FMR.concatFold $ FMR.mapReduceFold
    (MR.filterUnpack filterData)
    (FMR.assignKeysAndData @k)
    (FMR.foldAndAddKey $ weightedBinomialFold testData weightData)
-}

weightedCountFoldGeneral
  :: forall k r d
   . (Ord (F.Record k)
     , FI.RecVec (k V.++ CountCols)
     , k F.⊆ (k V.++ r)
     , d F.⊆ (k V.++ r)
     )
  => (F.Record r -> F.Record k)
  -> (F.Record r -> Bool) -- ^ include this row?
  -> (F.Record d -> Bool) -- ^ success ?
  -> (F.Record d -> Double) -- ^ weight
  -> FL.Fold (F.Record r) (F.FrameRec (k V.++ CountCols))
weightedCountFoldGeneral getKey filterData testData weightData =
  FL.prefilter filterData $ FMR.concatFold $ FMR.mapReduceFold
    (MR.Unpack $ \r ->  [getKey r `V.rappend` r])
    (FMR.assignKeysAndData @k @d)
    (FMR.foldAndAddKey $ weightedBinomialFold testData weightData)

-}

getFraction r =
  let n = F.rgetField @Count r
  in  if n == 0
        then 0
        else realToFrac (F.rgetField @Successes r) / realToFrac n

getFractionWeighted r =
  let n = F.rgetField @Count r
  in  if n == 0 then 0 else F.rgetField @WeightedSuccesses r / realToFrac n

data RecordColsProxy k = RecordColsProxy deriving (Show, Enum, Bounded, A.Ix, Eq, Ord)
type instance GLM.GroupKey (RecordColsProxy k) = F.Record k

--type instance GLM.GroupKey (Proxy k) = F.Record k

recordToGroupKey
  :: forall k r . (k F.⊆ r) => F.Record r -> RecordColsProxy k -> F.Record k
recordToGroupKey r _ = F.rcast @k r

glmErrorToPandocError :: GLM.GLMError -> PE.PandocError
glmErrorToPandocError x = PE.PandocSomeError $ show x

type GroupCols gs cs = gs V.++ cs
type CountedCols gs cs = GroupCols gs cs V.++ CountCols

---
newtype  SimplePredictor ps = SimplePredictor { unSimplePredictor :: F.Record ps }

instance (Show (F.Record ps)) => Show (SimplePredictor ps) where
  show (SimplePredictor x) = "SimplePredictor " ++ show @String x

instance (Eq (F.Record ps)) => Eq (SimplePredictor ps) where
  (SimplePredictor x) == (SimplePredictor y) = x == y

instance (Ord (F.Record ps)) => Ord (SimplePredictor ps) where
  compare (SimplePredictor x) (SimplePredictor y) = compare x y

instance (Ord (F.Record ps), K.FiniteSet (F.Record ps)) => Enum (SimplePredictor ps) where
  toEnum n =
    let im = IM.fromList $ zip [0..] $ SimplePredictor <$> S.toAscList K.elements
    in fromMaybe (error "Bad n given to toEnum :: Int -> SimplePredictor") $ IM.lookup n im
  fromEnum a =
    let m = M.fromList $ zip (SimplePredictor <$> S.toAscList K.elements) [0..]
    in fromMaybe (error "Bad SimplePredictor given to fromEnum :: SimplePredictor -> Int.  This means something is very bad.") $ M.lookup a m

instance (K.FiniteSet (F.Record ps)) => Bounded (SimplePredictor ps) where
  minBound = fromMaybe (error "Empty Set of Records used as predictor") $ viaNonEmpty head $ SimplePredictor <$> S.toList K.elements
  maxBound = fromMaybe (error "Empty Set of Records used as predictor") $ viaNonEmpty last $ SimplePredictor <$> S.toList K.elements

allSimplePredictors :: K.FiniteSet (F.Record ps) => [SimplePredictor ps]
allSimplePredictors = SimplePredictor <$> S.toList K.elements

type SimpleEffect ps = GLM.WithIntercept (SimplePredictor ps)

simplePredictor :: forall ps rs. (ps F.⊆ rs
                                     , Eq (F.Record ps)
                                     )
                    => F.Record rs -> SimplePredictor ps -> Double
simplePredictor r p = if F.rcast @ps r == unSimplePredictor p then 1 else 0

predMap :: forall cs. (K.FiniteSet (F.Record cs), Ord (F.Record cs), cs F.⊆ cs)
  => F.Record cs -> M.Map (SimplePredictor cs) Double
predMap r =  M.fromList $ Relude.fmapToSnd (simplePredictor r) allSimplePredictors

catPredMaps :: forall cs.  (K.FiniteSet (F.Record cs), Ord (F.Record cs), cs F.⊆ cs)
  => M.Map (F.Record cs) (M.Map (SimplePredictor cs) Double)
catPredMaps = M.fromList $ fmap (\k -> (unSimplePredictor k,predMap (unSimplePredictor k))) allSimplePredictors

data  LocationHolder c f a =  LocationHolder { locName :: T.Text
                                             , locKey :: Maybe (F.Rec f LocationCols)
                                             , catData :: M.Map (F.Rec f c) a
                                             } deriving (Generic)

deriving instance (V.RMap c
                  , V.ReifyConstraint Show F.ElField c
                  , V.RecordToList c
                  , Show a) => Show (LocationHolder c F.ElField a)

instance (SE.Serialize a
         , Ord (F.Rec FS.SElField c)
         , SE.GSerializePut
           (Rep (F.Rec FS.SElField c))
         , SE.GSerializeGet (Rep (F.Rec FS.SElField c))
         , (Generic (F.Rec FS.SElField c))
         ) => SE.Serialize (LocationHolder c FS.SElField a)

lhToS :: (Ord (F.Rec FS.SElField c)
         , V.RMap c
         )
      => LocationHolder c F.ElField a -> LocationHolder c FS.SElField a
lhToS (LocationHolder n lkM cdm) = LocationHolder n (fmap FS.toS lkM) (M.mapKeys FS.toS cdm)

lhFromS :: (Ord (F.Rec F.ElField c)
           , V.RMap c
         ) => LocationHolder c FS.SElField a -> LocationHolder c F.ElField a
lhFromS (LocationHolder n lkM cdm) = LocationHolder n (fmap FS.fromS lkM) (M.mapKeys FS.fromS cdm)

type LocationCols = '[BR.StateAbbreviation]
locKeyPretty :: F.Record LocationCols -> T.Text
locKeyPretty r =
  let stateAbbr = F.rgetField @BR.StateAbbreviation r
  in stateAbbr

--type ASER = '[BR.SimpleAgeC, BR.SexC, BR.CollegeGradC, BR.SimpleRaceC]
predictionsByLocation ::
  forall ps r rs. ( V.RMap ps
                  , V.ReifyConstraint Show V.ElField ps
                  , V.RecordToList ps
                  , Ord (F.Record ps)
                  , Ord (SimplePredictor ps)
                  , FI.RecVec (ps V.++ CountCols)
                  , F.ElemOf (ps V.++ CountCols) Count
                  , F.ElemOf (ps V.++ CountCols) MeanWeight
                  , F.ElemOf (ps V.++ CountCols) UnweightedSuccesses
                  , F.ElemOf (ps V.++ CountCols) VarWeight
                  , F.ElemOf (ps V.++ CountCols) WeightedSuccesses
                  , K.FiniteSet (F.Record ps)
                  , Show (F.Record (LocationCols V.++ ps V.++ CountCols))
                  , Show (F.Record ps)
                  , Enum (SimplePredictor ps)
                  , Bounded (SimplePredictor ps)
                  , (ps V.++ CountCols) F.⊆ (LocationCols V.++ ps V.++ CountCols)
                  , ps F.⊆ (ps V.++ CountCols)
                  , ps F.⊆ (LocationCols V.++ ps V.++ CountCols)
                  , F.ElemOf rs BR.StateAbbreviation
                  , K.KnitEffects r
               )
  => GLM.MinimizeDevianceVerbosity
  -> F.FrameRec rs
  -> FL.Fold (F.Record rs) (F.FrameRec (LocationCols V.++ ps V.++ CountCols))
  -> [SimpleEffect ps]
  -> M.Map (F.Record ps) (M.Map (SimplePredictor ps) Double)
  -> K.Sem r [LocationHolder ps V.ElField Double]
predictionsByLocation verbosity ccesFrame countFold predictors catPredMap =
  P.mapError glmErrorToPandocError  $ K.wrapPrefix "predictionsByLocation" $ do
    K.logLE K.Diagnostic "Inferring"
    (mm, rc, ebg, bu, vb, bs) <- inferMR @LocationCols @ps @ps
                                 verbosity
                                 countFold
                                 predictors
                                 simplePredictor
                                 ccesFrame

    let states = FL.fold FL.set $ fmap (F.rgetField @BR.StateAbbreviation) ccesFrame
        allStateKeys = (F.&: V.RNil) <$> FL.fold FL.list states
        predictLoc l = LocationHolder (locKeyPretty l) (Just l) catPredMap
        toPredict = [LocationHolder "National" Nothing catPredMap] <> fmap predictLoc allStateKeys
        predict (LocationHolder n lkM cpms) = P.mapError glmErrorToPandocError $ do
          let predictFrom catKey predMap =
                let groupKeyM = fmap (`V.rappend` catKey) lkM --lkM >>= \lk -> return $ lk `V.rappend` catKey
                    emptyAsNationalGKM = case groupKeyM of
                                           Nothing -> Nothing
                                           Just k -> k <$ GLM.categoryNumberFromKey rc k (RecordColsProxy @(LocationCols V.++ ps))
                in (fmap snd . GLM.runLogOnGLMException Nothing)
                   $ GLM.predictFromBetaB mm (`M.lookup` predMap) (const emptyAsNationalGKM) rc ebg bu vb
          cpreds <- M.traverseWithKey predictFrom cpms
          return $ LocationHolder n lkM cpreds
    traverse predict toPredict

---

-- TODO: Add bootstraps back in, make optional
-- Maybe set up way to serialize all this at this level??

-- ls: are keys for locations (groups with no matching fixed-effect)
-- cs: are keys for fixedEffects, one for each constructor in b
-- b: Fixed Effect type
-- g: group type
inferMR
  :: forall ls cs ks b g f rs effs
   . ( Foldable f
     , Functor f
     , K.KnitEffects effs
     , F.RDeleteAll ls (ls V.++ ks V.++ CountCols) ~ (ks V.++ CountCols)
     , (ls V.++ ks V.++ CountCols) ~ (ls V.++ (ks V.++ CountCols))
     , Ord (F.Record ls)
     , ls F.⊆ (ls V.++ ks V.++ CountCols)
     , (ks V.++ CountCols) F.⊆ (ls V.++ ks V.++ CountCols)
     , ks F.⊆ (ks V.++ CountCols)
     , (ls V.++ cs) F.⊆ (ls V.++ (ks V.++ CountCols))
     , Show (F.Record (ls V.++ ks V.++ CountCols))
     , FI.RecVec (ls V.++ (ks V.++ CountCols))
     , F.ElemOf (ks V.++ CountCols) Count
     , F.ElemOf (ks V.++ CountCols) MeanWeight
     , F.ElemOf (ks V.++ CountCols) VarWeight
     , F.ElemOf (ks V.++ CountCols) WeightedSuccesses
     , F.ElemOf (ks V.++ CountCols) UnweightedSuccesses
     , F.ElemOf (ls V.++ ks V.++ CountCols) Count
     , F.ElemOf (ls V.++ ks V.++ CountCols) MeanWeight
     , F.ElemOf (ls V.++ ks V.++ CountCols) VarWeight
     , F.ElemOf (ls V.++ ks V.++ CountCols) WeightedSuccesses
     , F.ElemOf (ls V.++ ks V.++ CountCols) UnweightedSuccesses
     , Show (F.Record (ls V.++ cs))
     , K.FiniteSet (F.Record ks)
     , Ord b
     , Show b
     , Enum b
     , Bounded b
     , Ord (F.Record (GroupCols ls cs))
     , g ~ RecordColsProxy (GroupCols ls cs)
     )
  => GLM.MinimizeDevianceVerbosity
  -> FL.Fold (F.Record rs) (F.FrameRec (ls V.++ ks V.++ CountCols))
  -> [GLM.WithIntercept b] -- fixed effects to fit
  -> (F.Record (ls V.++ ks V.++ CountCols) -> b -> Double)  -- how to get a fixed effect value from a record
  -> f (F.Record rs)
  -> K.Sem
       effs
       ( GLM.MixedModel b g
       , GLM.RowClassifier g
       , GLM.EffectsByGroup g b
       , GLM.BetaVec
       , LA.Vector Double
       , [(GLM.BetaVec, LA.Vector Double)]
       )
inferMR verbosity cf fixedEffectList getFixedEffect rows =
  P.mapError glmErrorToPandocError
  $ (fmap snd . GLM.runLogOnGLMException Nothing)
  $ K.wrapPrefix "inferMR"
  $ do
    let
      addZeroCountsF = FMR.concatFold $ FMR.mapReduceFold
                       FMR.noUnpack
                       (FMR.splitOnKeys @ls)
                       ( FMR.makeRecsWithKey id
                         $ FMR.ReduceFold
                         $ const
                         $ K.addDefaultRec @ks zeroCount
                       )
      counted = FL.fold FL.list $ FL.fold addZeroCountsF $ FL.fold cf rows
--    K.logLE K.Diagnostic $ T.intercalate "\n" $ fmap (T.pack . show) counted
    let
      vCounts = VS.fromList $ fmap (F.rgetField @Count) counted
      designEffect mw vw = 1 + (vw / (mw * mw))
      vWeights = VS.fromList $ fmap
        (\r ->
            let mw = F.rgetField @MeanWeight r
                vw = F.rgetField @VarWeight r
            in  case (mw < 1e-12, vw < 1e-12) of
              (True, _) -> 0
              (False, True) -> 1
              _ -> 1/sqrt (designEffect mw vw)
        )
        counted -- VS.replicate (VS.length vCounts) 1.0

      fixedEffects = GLM.FixedEffects $ IS.fromList fixedEffectList
      groups       = IS.fromList [RecordColsProxy]
      (observations, fixedEffectsModelMatrix, rcM) = FL.fold
        (lmePrepFrame getFractionWeighted
         fixedEffects
         groups
         getFixedEffect
         (recordToGroupKey @(GroupCols ls cs))
        )
        counted
--        K.logLE K.Diagnostic $ "vCounts=" <> (T.pack $ show vCounts)
--        K.logLE K.Diagnostic $ "vWeights=" <> (T.pack $ show vWeights)
--        K.logLE K.Diagnostic $ "mX=" <> (T.pack $ show fixedEffectsModelMatrix)
--        K.logLE K.Diagnostic $ "vY=" <> (T.pack $ show observations)
    let regressionModelSpec = GLM.RegressionModelSpec
          fixedEffects
          fixedEffectsModelMatrix
          observations
    rowClassifier <- case rcM of
      Left  msg -> K.knitError msg
      Right x   -> return x
--        K.logLE K.Diagnostic $ "rc=" <> (T.pack $ show rowClassifier)
    let effectsByGroup =
          M.fromList [(RecordColsProxy, IS.fromList [GLM.Intercept])]
    fitSpecByGroup <- GLM.fitSpecByGroup @b @g fixedEffects
                      effectsByGroup
                      rowClassifier
    let lmmControls = GLM.LMMControls GLM.LMM_BOBYQA 1e-6 Nothing
        lmmSpec     = GLM.LinearMixedModelSpec
                      (GLM.MixedModelSpec regressionModelSpec fitSpecByGroup)
                      lmmControls
        cc = GLM.PIRLSConvergenceCriterion GLM.PCT_Deviance 1e-6 20
        glmmControls = GLM.GLMMControls GLM.UseCanonical 10 cc
        glmmSpec = GLM.GeneralizedLinearMixedModelSpec
                   lmmSpec
                   vWeights
                   (GLM.Binomial vCounts)
                   glmmControls
        mixedModel = GLM.GeneralizedLinearMixedModel glmmSpec
    randomEffectsModelMatrix <- GLM.makeZ fixedEffectsModelMatrix
                                fitSpecByGroup
                                rowClassifier
--        K.logLE K.Diagnostic $ "smZ=" <> (T.pack $ show randomEffectsModelMatrix)
    let randomEffectCalc = GLM.RandomEffectCalculated
          randomEffectsModelMatrix
          (GLM.makeLambda fitSpecByGroup)
        th0         = GLM.setCovarianceVector fitSpecByGroup 1 0
        -- mdVerbosity = MDVNone
    GLM.checkProblem mixedModel randomEffectCalc
    K.logLE K.Info "Fitting data..."
    ((th, pd, sigma2, beta, mzvu, mzvb, cs), vMuSol, cf) <- GLM.minimizeDeviance
                                                            verbosity
                                                            ML
                                                            mixedModel
                                                            randomEffectCalc
                                                            th0
    vb <- K.knitMaybe "b-vector (random effects coefficients) is zero in MRP.inferMR" $ GLM.useMaybeZeroVec Nothing Just mzvb
    GLM.report mixedModel
      randomEffectsModelMatrix
      beta
      (SD.toSparseVector vb)
    let fes = GLM.fixedEffectStatistics mixedModel sigma2 cs beta
    K.logLE K.Diagnostic $ "FixedEffectStatistics: " <> show fes
    epg <- GLM.effectParametersByGroup @g @b rowClassifier effectsByGroup vb
--        K.logLE K.Diagnostic
--          $  "EffectParametersByGroup: "
--          <> (T.pack $ show epg)
    gec <- GLM.effectCovariancesByGroup effectsByGroup mixedModel sigma2 th
    K.logLE K.Diagnostic
      $  "EffectCovariancesByGroup: "
      <> show gec
    rebl <- GLM.randomEffectsByLabel epg rowClassifier
    K.logLE K.Diagnostic
      $  "Random Effects:\n"
      <> GLM.printRandomEffectsByLabel rebl
    smCondVar <- GLM.conditionalCovariances mixedModel
                 cf
                 randomEffectCalc
                 th
                 beta
                 mzvu
    let bootstraps                           = []
    let GLM.FixedEffectStatistics _ mBetaCov = fes
    let f r = do
                let obs = getFractionWeighted r
                predictCVCI <- GLM.predictWithCI
                  mixedModel
                  (Just . getFixedEffect r)
                  (Just . recordToGroupKey @(GroupCols ls cs) r)
                  rowClassifier
                  effectsByGroup
                  beta
                  vb
                  (ST.mkCL 0.95)
                  (GLM.NaiveCondVarCI mBetaCov smCondVar)
                return (r, obs, predictCVCI)
    fitted <- traverse f (FL.fold FL.list counted)
    K.logLE K.Diagnostic
      $  "Fitted:\n"
      <> (T.intercalate "\n" $ fmap show fitted)
    fixedEffectTable <- GLM.printFixedEffects fes
    K.logLE K.Diagnostic $ "FixedEffects:\n" <> fixedEffectTable
    let GLM.FixedEffectStatistics fep _ = fes
    return
      (mixedModel, rowClassifier, effectsByGroup, beta, vb, bootstraps) -- fes, epg, rowClassifier, bootstraps)

lmePrepFrame
  :: forall p g rs
   . (GLM.PredictorC p, GLM.GroupC g)
  => (F.Record rs -> Double) -- ^ observations
  -> GLM.FixedEffects p
  -> IS.IndexedSet g
  -> (F.Record rs -> p -> Double) -- ^ predictors
  -> (F.Record rs -> g -> GLM.GroupKey g)  -- ^ classifiers
  -> FL.Fold
       (F.Record rs)
       ( LA.Vector Double
       , LA.Matrix Double
       , Either T.Text (GLM.RowClassifier g)
       ) -- ^ (X,y,(row-classifier, size of class))
lmePrepFrame observationF fe groupIndices getPredictorF classifierLabelF
  = let
      makeInfoVector
        :: M.Map g (M.Map (GLM.GroupKey g) Int)
        -> M.Map g (GLM.GroupKey g)
        -> Either T.Text (V.Vector (GLM.ItemInfo g))
      makeInfoVector indexMaps groupKeys =
        let
          g (grp, groupKey) =
            GLM.ItemInfo
              <$> (   maybe
                      (Left $ "Failed on " <> show (grp, groupKey))
                      Right
                  $   M.lookup grp indexMaps
                  >>= M.lookup groupKey
                  )
              <*> pure groupKey
        in  fmap V.fromList $ traverse g $ M.toList groupKeys
      makeRowClassifier
        :: Traversable f
        => M.Map g (M.Map (GLM.GroupKey g) Int)
        -> f (M.Map g (GLM.GroupKey g))
        -> Either T.Text (GLM.RowClassifier g)
      makeRowClassifier indexMaps labels = do
        let sizes = fmap M.size indexMaps
        indexed <- traverse (makeInfoVector indexMaps) labels
        return $ GLM.RowClassifier groupIndices
                                   sizes
                                   (V.fromList $ FL.fold FL.list indexed)
                                   indexMaps
      getPredictorF' _   GLM.Intercept     = 1
      getPredictorF' row (GLM.Predictor x) = getPredictorF row x
      predictorF row = LA.fromList $ case fe of
        GLM.FixedEffects indexedFixedEffects ->
          getPredictorF' row <$> IS.members indexedFixedEffects
        GLM.InterceptOnly -> [1]
      getClassifierLabels :: F.Record rs -> M.Map g (GLM.GroupKey g)
      getClassifierLabels r =
        M.fromList $ Relude.fmapToSnd (classifierLabelF r) $ IS.members groupIndices
      foldObs   = LA.fromList <$> FL.premap observationF FL.list
      foldPred  = LA.fromRows <$> FL.premap predictorF FL.list
      foldClass = FL.premap getClassifierLabels FL.list
      g (vY, mX, ls) =
        ( vY
        , mX
        , makeRowClassifier
          (snd $ State.execState (addAll ls) (M.empty, M.empty))
          ls
        )
    in
      g <$> ((,,) <$> foldObs <*> foldPred <*> foldClass)

addOne
  :: GLM.GroupC g
  => (g, GLM.GroupKey g)
  -> State.State (M.Map g Int, M.Map g (M.Map (GLM.GroupKey g) Int)) ()
addOne (grp, label) = do
  (nextIndexMap, groupIndexMaps) <- State.get
  let groupIndexMap = fromMaybe M.empty $ M.lookup grp groupIndexMaps
  case M.lookup label groupIndexMap of
    Nothing -> do
      let index         = fromMaybe 0 $ M.lookup grp nextIndexMap
          nextIndexMap' = M.insert grp (index + 1) nextIndexMap
          groupIndexMaps' =
            M.insert grp (M.insert label index groupIndexMap) groupIndexMaps
      State.put (nextIndexMap', groupIndexMaps')
      return ()
    _ -> return ()

addMany
  :: (GLM.GroupC g, Traversable h)
  => h (g, GLM.GroupKey g)
  -> State.State (M.Map g Int, M.Map g (M.Map (GLM.GroupKey g) Int)) ()
addMany = traverse_ addOne

addAll
  :: GLM.GroupC g
  => [M.Map g (GLM.GroupKey g)]
  -> State.State (M.Map g Int, M.Map g (M.Map (GLM.GroupKey g) Int)) ()
addAll = traverse_ (addMany . M.toList)

-- useful data folds
weightedSumF :: (Real w, Real v, Fractional y) => FL.Fold (w, v) y
weightedSumF =
  (/)
    <$> FL.premap (\(w, v) -> realToFrac w * realToFrac v) FL.sum
    <*> P.dimap fst realToFrac FL.sum

weightedSumRecF
  :: forall w v y rs
   . ( V.KnownField w
     , V.KnownField v
     , Fractional y
     , F.ElemOf rs w
     , F.ElemOf rs v
     , Real (V.Snd w)
     , Real (V.Snd v)
     )
  => FL.Fold (F.Record rs) y
weightedSumRecF =
  FL.premap (\r -> (F.rgetField @w r, F.rgetField @v r)) weightedSumF

sumProdIfRecF
  :: forall w v t y rs
   . ( V.KnownField w
     , V.KnownField v
     , V.KnownField t
     , Fractional y
     , F.ElemOf rs w
     , F.ElemOf rs v
     , F.ElemOf rs t
     , Real (V.Snd w)
     , Real (V.Snd v)
     )
  => (V.Snd t -> Bool)
  -> FL.Fold (F.Record rs) y
sumProdIfRecF test = FL.prefilter (test . F.rgetField @t) $ FL.premap
  (\r -> realToFrac (F.rgetField @w r) * realToFrac (F.rgetField @v r))
  FL.sum

max_stars_repo_name: blueripple/preference-model

max_issues_repo_name: blueripple/preference-model

max_forks_repo_name: blueripple/preference-model


text:
module Main where

import App
import UKF

import Numeric.LinearAlgebra ((!))

app conn = do
  sensor <- recv conn
  app' conn (makeFilter sensor)

app' conn kf = do
  sensor <- recv conn

  let kf' = update sensor kf
  let px = (kf_x kf') ! 0
  let py = (kf_x kf') ! 1

  send conn Response { rX = px
                     , rY = py }
  app' conn kf'

main =
  runApp app

ext: hs


text:
{-# LANGUAGE CPP                                      #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TemplateHaskell                          #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeOperators                            #-}
{-# OPTIONS_GHC -fno-warn-orphans                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

#if MIN_VERSION_base(4,12,0)
{-# LANGUAGE NoStarIsType         #-}
#endif

import           Control.Monad
import           Data.Complex
import           Data.Function
import           Data.Functor.Identity
import           Data.Maybe
import           Data.Proxy
import           GHC.TypeLits
import           Hedgehog
import           Numeric.Natural
import           System.Exit
import           System.IO
import qualified Data.Vector.Generic                 as UVG
import qualified Data.Vector.Generic.Sized           as VG
import qualified Data.Vector.Sized                   as V
import qualified Data.Vector.Storable.Sized          as VS
import qualified Hedgehog.Gen                        as Gen
import qualified Hedgehog.Range                      as Range
import qualified Numeric.LinearAlgebra               as HU
import qualified Numeric.LinearAlgebra.Static        as H
import qualified Numeric.LinearAlgebra.Static.Vector as H

instance KnownNat n => Eq (H.R n) where
    (==) = (==) `on` H.extract

instance KnownNat n => Eq (H.C n) where
    (==) = (==) `on` H.extract

instance (KnownNat m, KnownNat n) => Eq (H.L m n) where
    (==) = (==) `on` H.extract

instance (KnownNat m, KnownNat n) => Eq (H.M m n) where
    (==) = (==) `on` H.extract

someNatVal' :: Natural -> SomeNat
someNatVal' = fromJust . someNatVal . fromIntegral

genDouble :: Gen H.ℝ
genDouble = Gen.double (Range.linearFracFrom 0 (-10) 10)

genComplex :: Gen H.ℂ
genComplex = (:+) <$> genDouble <*> genDouble

sizeRange :: Num a => Range a
sizeRange = Range.constant 5 10

validVector
    :: (UVG.Vector v a, KnownNat n, MonadTest m)
    => VG.Vector v n a
    -> m ()
validVector v = VG.length v === UVG.length (VG.fromSized v)

prop_rVec :: Property
prop_rVec = property $ do
    xs <- forAll $ Gen.list sizeRange genDouble
    SomeNat (Proxy :: Proxy n) <- pure $ someNatVal' (fromIntegral (length xs))
    tripping (H.vector @n xs)
             (VS.map (* 2) . H.rVec)
             (Identity . (/ 2) . H.vecR)

prop_rVec_konst :: Property
prop_rVec_konst = property $ do
    x <- forAll genDouble
    SomeNat (Proxy :: Proxy n) <- fmap someNatVal' . forAll
                                $ Gen.integral sizeRange
    validVector . H.rVec @n $ H.konst x

prop_vecR :: Property
prop_vecR = property $ do
    xs <- forAll $ Gen.list sizeRange genDouble
    VS.withSizedList xs $ \v ->
      tripping v ((* 2) . H.vecR)
                 (Identity . VS.map (/ 2) . H.rVec)

prop_cVec :: Property
prop_cVec = property $ do
    xs <- forAll $ Gen.list sizeRange genComplex
    SomeNat (Proxy :: Proxy n) <- pure $ someNatVal' (fromIntegral (length xs))
    tripping (fromJust . H.create $ HU.fromList xs)
             (VS.map (* 2) . H.cVec @n)
             (Identity . (/ 2) . H.vecC)

prop_cVec_konst :: Property
prop_cVec_konst = property $ do
    x <- forAll genComplex
    SomeNat (Proxy :: Proxy n) <- fmap someNatVal' . forAll
                                $ Gen.integral sizeRange
    validVector . H.cVec @n $ H.konst x

prop_vecC :: Property
prop_vecC = property $ do
    xs <- forAll $ Gen.list sizeRange genComplex
    VS.withSizedList xs $ \v ->
      tripping v ((* 2) . H.vecC)
                 (Identity . VS.map (/ 2) . H.cVec)

genMatList :: Gen a -> Gen (SomeNat, SomeNat, [[a]])
genMatList g = do
    m  <- Gen.integral sizeRange
    n  <- Gen.integral sizeRange
    xs <- (replicateM (fromIntegral m) . replicateM (fromIntegral n)) g
    return ( someNatVal' m
           , someNatVal' n
           , xs
           )

genMatKonst :: Gen a -> Gen (SomeNat, SomeNat, a)
genMatKonst g = do
    m <- Gen.integral sizeRange
    n <- Gen.integral sizeRange
    x <- g
    return ( someNatVal' m
           , someNatVal' n
           , x
           )

prop_lRows :: Property
prop_lRows = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , xs
     ) <- forAll $ genMatList genDouble
    tripping (fromJust . H.create $ HU.fromLists xs)
             (V.map (* 2) . H.lRows @m @n)
             (Identity . (/ 2) . H.rowsL)

prop_lRows_konst :: Property
prop_lRows_konst = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , x
     ) <- forAll $ genMatKonst genDouble
    validVector . H.lRows @m @n $ H.konst x

prop_rowsL :: Property
prop_rowsL = property $ do
    (_, SomeNat (Proxy :: Proxy n), xs)  <- forAll $ genMatList genDouble
    V.withSizedList xs $ \(v :: V.Vector m [Double]) ->
      tripping (V.map (fromJust . H.create . HU.fromList) v)
               ((* 2) . H.rowsL @m @n)
               (Identity . V.map (/ 2) . H.lRows)

prop_rowsL_konst :: Property
prop_rowsL_konst = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , x
     ) <- forAll $ genMatKonst genDouble
    let rs = V.replicate $ H.konst x
    validVector . H.lVec . H.rowsL @m @n $ rs
    tripping rs
        ((* 2) . H.rowsL @m @n)
        (Identity . V.map (/ 2) . H.lRows)

prop_lCols :: Property
prop_lCols = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , xs
     ) <- forAll $ genMatList genDouble
    tripping (fromJust . H.create $ HU.fromLists xs)
             (V.map (* 2) . H.lCols @m @n)
             (Identity . (/ 2) . H.colsL)

prop_lCols_konst :: Property
prop_lCols_konst = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , x
     ) <- forAll $ genMatKonst genDouble
    validVector . H.lCols @m @n $ H.konst x

prop_colsL :: Property
prop_colsL = property $ do
    (_, SomeNat (Proxy :: Proxy m), xs)  <- forAll $ genMatList genDouble
    V.withSizedList xs $ \(v :: V.Vector n [Double]) ->
      tripping (V.map (fromJust . H.create . HU.fromList) v)
               ((* 2) . H.colsL @m @n)
               (Identity . V.map (/ 2) . H.lCols)

prop_colsL_konst :: Property
prop_colsL_konst = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , x
     ) <- forAll $ genMatKonst genDouble
    let cs = V.replicate $ H.konst x
    validVector . H.lVec . H.colsL @m @n $ cs
    tripping cs
        ((* 2) . H.colsL @m @n)
        (Identity . V.map (/ 2) . H.lCols)

prop_vecL :: Property
prop_vecL = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , xs
     ) <- forAll $ genMatList genDouble
    let v :: VS.Vector (m * n) Double
        v = fromJust $ VS.fromList @_ @(m * n) (concat xs)
    tripping v ((* 2) . H.vecL @m @n)
              (Identity . VS.map (/ 2) . H.lVec)

prop_lVec :: Property
prop_lVec = property $ do
    ( SomeNat (pM@Proxy :: Proxy m)
     , SomeNat (pN@Proxy :: Proxy n)
     , xs
     ) <- forAll $ genMatList genDouble
    let m :: H.L m n
        m = fromJust
          . H.create
          . (fromIntegral (natVal pM) HU.>< fromIntegral (natVal pN))
          $ concat xs
    tripping m (VS.map (* 2) . H.lVec)
              (Identity . (/ 2) . H.vecL)

prop_lVec_konst :: Property
prop_lVec_konst = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , x
     ) <- forAll $ genMatKonst genDouble
    let xs = H.konst x
    validVector $ H.lVec @m @n xs
    tripping xs
        (VS.map (* 2) . H.lVec @m @n)
        (Identity . (/ 2) . H.vecL)

prop_mRows :: Property
prop_mRows = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , xs
     ) <- forAll $ genMatList genComplex
    tripping (fromJust . H.create $ HU.fromLists xs)
             (V.map (* 2) . H.mRows @m @n)
             (Identity . (/ 2) . H.rowsM)

prop_mRows_konst :: Property
prop_mRows_konst = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , x
     ) <- forAll $ genMatKonst genComplex
    validVector . H.mRows @m @n $ H.konst x

prop_rowsM :: Property
prop_rowsM = property $ do
    (_, SomeNat (Proxy :: Proxy n), xs)  <- forAll $ genMatList genComplex
    V.withSizedList xs $ \(v :: V.Vector m [Complex Double]) ->
      tripping (V.map (fromJust . H.create . HU.fromList) v)
               ((* 2) . H.rowsM @m @n)
               (Identity . V.map (/ 2) . H.mRows)

prop_rowsM_konst :: Property
prop_rowsM_konst = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , x
     ) <- forAll $ genMatKonst genComplex
    let rs = V.replicate $ H.konst x
    validVector . H.mVec . H.rowsM @m @n $ rs
    tripping rs
        ((* 2) . H.rowsM @m @n)
        (Identity . V.map (/ 2) . H.mRows)

prop_mCols :: Property
prop_mCols = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , xs
     ) <- forAll $ genMatList genComplex
    tripping (fromJust . H.create $ HU.fromLists xs)
             (V.map (* 2) . H.mCols @m @n)
             (Identity . (/ 2) . H.colsM)

prop_mCols_konst :: Property
prop_mCols_konst = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , x
     ) <- forAll $ genMatKonst genComplex
    validVector . H.mCols @m @n $ H.konst x

prop_colsM :: Property
prop_colsM = property $ do
    (_, SomeNat (Proxy :: Proxy m), xs)  <- forAll $ genMatList genComplex
    V.withSizedList xs $ \(v :: V.Vector n [Complex Double]) ->
      tripping (V.map (fromJust . H.create . HU.fromList) v)
               ((* 2) . H.colsM @m @n)
               (Identity . V.map (/ 2) . H.mCols)

prop_colsM_konst :: Property
prop_colsM_konst = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , x
     ) <- forAll $ genMatKonst genComplex
    let cs = V.replicate $ H.konst x
    validVector . H.mVec . H.colsM @m @n $ cs
    tripping cs
        ((* 2) . H.colsM @m @n)
        (Identity . V.map (/ 2) . H.mCols)

prop_vecM :: Property
prop_vecM = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , xs
     ) <- forAll $ genMatList genComplex
    let v :: VS.Vector (m * n) (Complex Double)
        v = fromJust $ VS.fromList @_ @(m * n) (concat xs)
    tripping v ((* 2) . H.vecM @m @n)
              (Identity . VS.map (/ 2) . H.mVec)

prop_mVec :: Property
prop_mVec = property $ do
    ( SomeNat (pM@Proxy :: Proxy m)
     , SomeNat (pN@Proxy :: Proxy n)
     , xs
     ) <- forAll $ genMatList genComplex
    let m :: H.M m n
        m = fromJust
          . H.create
          . (fromIntegral (natVal pM) HU.>< fromIntegral (natVal pN))
          $ concat xs
    tripping m (VS.map (* 2) . H.mVec)
               (Identity . (/ 2) . H.vecM)

prop_mVec_konst :: Property
prop_mVec_konst = property $ do
    ( SomeNat (Proxy :: Proxy m)
     , SomeNat (Proxy :: Proxy n)
     , x
     ) <- forAll $ genMatKonst genComplex
    let xs = H.konst x
    validVector $ H.mVec @m @n xs
    tripping xs
        (VS.map (* 2) . H.mVec @m @n)
        (Identity . (/ 2) . H.vecM)

main :: IO ()
main = do
  hSetBuffering stdout LineBuffering
  hSetBuffering stderr LineBuffering

  results <- checkParallel $$(discover)

  unless results exitFailure

max_stars_repo_name: mstksg/hmatrix-vector-sized

max_issues_repo_name: mstksg/hmatrix-vector-sized

max_forks_repo_name: mstksg/hmatrix-vector-sized


text:
-- TODO: better naming conventions

{-# LANGUAGE TypeFamilies, TypeOperators #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE DataKinds, PolyKinds, ConstraintKinds #-}
{-# LANGUAGE FlexibleInstances, FlexibleContexts, UndecidableInstances #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PartialTypeSignatures #-}

module GrammarNNDAG where

import Control.Monad (foldM)
import Data.Function (fix)

import Generics.SOP hiding (SingI)
import Generics.SOP.NP
import Generics.SOP.NS
import Generics.SOP.Constraint
import Generics.SOP.Dict
import Data.Vinyl hiding (Dict)
import Data.Singletons.Prelude hiding (All, And)
import Data.Singletons.Prelude.List (Length)

import Utils
import List
import Grammar
import TensorHMatrix
import DAGIO
import TensorDAG
import VinylUtils
import Random

import Data.Default
import DefaultM

import Numeric.LinearAlgebra (Numeric)
import qualified Data.Vector.Storable as Vector

type family Size (s :: k -> *) t :: k
type IntegralS (s :: k -> *) = IntegralK ('KProxy :: KProxy k)

class SingI (Size s t) => KnownSize s t
instance SingI (Size s t) => KnownSize s t

class All2 (KnownSize s) (Code t) => KnownSizes s t
instance All2 (KnownSize s) (Code t) => KnownSizes s t

class SingI (FromInteger (Length (Code t)) :: k) => KnownCode (s :: k -> *) t
instance SingI (FromInteger (Length (Code t)) :: k) => KnownCode (s :: k -> *) t

class HasNodes t where
  getNodes :: t -> [Some Node]

newtype Repr a dim = Repr { runRepr :: Node (Tensor '[dim] a) }

instance HasNodes (Repr a dim) where
  getNodes (Repr node) = [Some node]

newtype ReprT a s t = ReprT { runReprT :: Node (Tensor '[Size s t] a) }

instance HasNodes (ReprT a s t) where
  getNodes (ReprT node) = [Some node]

instance (Default a, Usable a, IntegralS s, KnownSize s t) => DefaultM IO (ReprT a s t) where
  defM = ReprT <$> defM

--newtype Linear a outDim inDim = Linear (Node (Tensor '[outDim, inDim] a))
newtype LinearT a s tOut tIn = LinearT (Node (Tensor '[Size s tOut, Size s tIn] a))

instance HasNodes (LinearT a s tOut tIn) where
  getNodes (LinearT node) = [Some node]

--instance (Default a, Usable a, IntegralN outDim, SingI outDim, SingI inDim) => DefaultM IO (Linear a outDim inDim) where
--  defM = Linear <$> defM

instance (Default a, Usable a, IntegralS s, KnownSize s tIn, KnownSize s tOut) => DefaultM IO (LinearT a s tOut tIn) where
  defM = LinearT <$> defM

--linear :: Numeric a => Linear a outDim inDim -> Repr a inDim -> IO (Node (Tensor a '[outDim]))
--linear (Linear m) (Repr v) = makeMV m v

linearT :: (SingI (Size s tOut), IntegralS s, Usable a) => LinearT a s tOut tIn -> ReprT a s tIn -> IO (ReprT a s tOut)
linearT (LinearT m) (ReprT v) = ReprT <$> makeMV m v

data EncodeParams' a s parent children = EncodeParams' (ReprT a s parent) (NP (LinearT a s parent) children)

instance HasNodes (EncodeParams' a s parent children) where
  getNodes (EncodeParams' bias weights) = getNodes bias ++ concat (collapse_NP $ liftA_NP' (K . getNodes) weights)

instance (Default a, Usable a, IntegralS s, KnownSize s parent, SListI children, All (KnownSize s) children) => DefaultM IO (EncodeParams' a s parent children) where
  defM = EncodeParams' <$> defM <*> sequence'_NP (cpure_NP (Proxy::Proxy (KnownSize s)) (Comp defM))

encodeParent' :: (Usable a, Floating a, Floating (Vector.Vector a), IntegralS s, KnownSize s parent) =>
  EncodeParams' a s parent children -> NP (ReprT a s) children -> IO (ReprT a s parent)
encodeParent' (EncodeParams' (ReprT bias) weights) inputs = do
  let linearT' m v = Comp $ K <$> linearT m v
  children <- map runReprT . collapse_NP <$> (sequence'_NP $ liftA2_NP' linearT' weights inputs)
  parent <- foldM (makeBinary (+)) bias children
  ReprT <$> makeUnary tanh parent

newtype EncodeParams a s t = EncodeParams { runEncodeParams :: NP (EncodeParams' a s t) (Code t) }

instance HasNodes (EncodeParams a s t) where
  getNodes (EncodeParams params) = concat . collapse_NP $ liftA_NP' (K . getNodes) params

encodeParent :: (Floating a, Floating (Vector.Vector a), Usable a, IntegralS s, KnownSize s t) =>
  EncodeParams a s t -> SOP (ReprT a s) (Code t) -> IO (ReprT a s t)
encodeParent (EncodeParams params) (SOP sop) = collapse_NS $ liftA2_NS' encode params sop
  where encode p cs = K $ encodeParent' p cs

data Encoding a s t where
  Primitive :: ReprT a s t -> Encoding a s t
  Generic :: ReprT a s t -> SOP (Encoding a s) (Code t) -> Encoding a s t

getRepr :: Encoding a s t -> ReprT a s t
getRepr (Primitive repr) = repr
getRepr (Generic repr _) = repr

instance (Default a, Usable a, IntegralS s, KnownSize s t, KnownSizes s t) => DefaultM IO (EncodeParams a s t) where
  defM = EncodeParams <$> sequence'_NP (cpure_NP (Proxy::Proxy (All (KnownSize s))) (Comp defM))

newtype Encoder a s t = Encoder { runEncoder :: t -> IO (Encoding a s t) }

encodeRec :: forall ts a s t. (Floating a, Floating (Vector.Vector a), Usable a, KnownSize s t, IntegralS s) =>
  Rec (Encoder a s) ts -> Dict (Contained ts) t -> EncodeParams a s t -> Encoder a s t

encodeRec encoders Dict params = Encoder f where
  encoders' = cpure_POP (Proxy::Proxy (Find ts)) (Fn $ Comp . runEncoder (index find encoders) . unI)
  f t = do
    children <- sequence_SOP' (ap_SOP encoders' (from t))
    let childReprs = liftA_SOP' getRepr children
    parent <- encodeParent params childReprs
    return $ Generic parent children

makeEncoders :: forall p g a s. (Floating a, Floating (Vector.Vector a), Usable a, IntegralS s, All (KnownSize s) g) =>
  Complete p g -> NP (EncodeParams a s) g -> Rec (Encoder a s) p -> Rec (Encoder a s) (p :++ g)
makeEncoders complete params prim = fix f where
  f encoders = rAppend prim $ np2rec (cliftA2_NP (Proxy::Proxy (KnownSize s)) (encodeRec encoders) (unAll_NP complete) params)

newtype AnyEncoder a s ts = AnyEncoder { runAnyEncoder :: forall t. Find ts t => t -> IO (Encoding a s t) }

makeEncoder :: (Floating a, Floating (Vector.Vector a), Usable a, IntegralS s, All (KnownSize s) g) =>
  Complete p g -> NP (EncodeParams a s) g -> Rec (Encoder a s) p -> AnyEncoder a s (p :++ g)
makeEncoder complete params prim = AnyEncoder $ runEncoder (index find encoders)
  where encoders = makeEncoders complete params prim

newtype LinearIn a s t t' = LinearIn (Node (Tensor '[Size s t', Size s t] a))

instance HasNodes (LinearIn a s t t') where
  getNodes (LinearIn node) = [Some node]

instance (Default a, Usable a, IntegralS s, KnownSize s t, KnownSize s t') => DefaultM IO (LinearIn a s t t') where
  defM = LinearIn <$> defM

linearIn :: (IntegralS s, KnownSize s t', Usable a) => LinearIn a s t t' -> ReprT a s t -> IO (ReprT a s t')
linearIn (LinearIn m) (ReprT v) = ReprT <$> makeMV m v

data AffineIn a s t t' = AffineIn (ReprT a s t') (LinearIn a s t t')

instance HasNodes (AffineIn a s t t') where
  getNodes (AffineIn bias linear) = getNodes bias ++ getNodes linear

instance (Default a, Usable a, IntegralS s, KnownSize s t, KnownSize s t') => DefaultM IO (AffineIn a s t t') where
  defM = AffineIn <$> defM <*> defM

affineIn :: (Floating a, Floating (Vector.Vector a), Usable a, IntegralS s, KnownSize s t') => AffineIn a s t t' -> ReprT a s t -> IO (ReprT a s t')
affineIn (AffineIn (ReprT bias) weights) input = do
  ReprT l <- linearIn weights input
  b <- makeBinary (+) bias l
  t <- makeUnary tanh b
  return $ ReprT t

data Affine' a outDim inDim =
  Affine' (Node (Tensor '[outDim] a)) (Node (Tensor '[outDim, inDim] a))

instance HasNodes (Affine' a outDim inDim) where
  getNodes (Affine' bias linear) = [Some bias, Some linear]

instance (Default a, Usable a, IntegralN outDim, SingI outDim, SingI inDim) => DefaultM IO (Affine' a outDim inDim) where
  defM = Affine' <$> defM <*> defM

sigmoid x = 1 / (1 + exp (-x))

affine' :: (Floating a, Floating (Vector.Vector a), Usable a, IntegralN outDim, SingI outDim) => Affine' a outDim inDim -> Repr a inDim -> IO (Repr a outDim)
affine' (Affine' bias weight) (Repr input) = do
  l <- makeMV weight input
  b <- makeBinary (+) bias l
  Repr <$> makeUnary sigmoid b

data DecodeParams a s t =
  DecodeParams (Affine' a (Length (Code t)) (Size s t)) (POP (AffineIn a s t) (Code t))

instance HasNodes (DecodeParams a s t) where
  getNodes (DecodeParams aff params) = getNodes aff ++ concat (collapse_POP' $ liftA_POP' (K . getNodes) params)

instance (Default a, Usable a, KnownCode s t, KnownSize s t, KnownSizes s t) => DefaultM IO (DecodeParams a s t) where
  defM = DecodeParams <$> defM <*> sequence'_POP (cpure_POP (Proxy::Proxy (KnownSize s)) (Comp defM))

decodeParent :: forall a s t. (Real a, Floating a, Floating (Vector.Vector a), Usable a, All2 (KnownSize s) (Code t), KnownCode s t) =>
  DecodeParams a s t -> ReprT a s t -> IO (SOP (ReprT a s) (Code t))
decodeParent (DecodeParams aff params) parent = do
  Repr node <- affine' aff (Repr $ runReprT parent)
  Vector v <- evalNode node
  let weights = map toRational $ Vector.toList v

  let children = np2ns . unPOP $ cliftA_POP (Proxy::Proxy (KnownSize s)) (Comp . flip affineIn parent) params

  child <- sample $ zip children weights
  sequence'_SOP $ SOP child

newtype Decoder a s t = Decoder { runDecoder :: ReprT a s t -> IO t }

decodeRec :: forall ts a s t. (Real a, Floating a, Floating (Vector.Vector a), Usable a, KnownSizes s t, KnownCode s t) =>
  Rec (Decoder a s) ts -> Dict (Contained ts) t -> DecodeParams a s t -> Decoder a s t

decodeRec decoders Dict params = Decoder f where
  decoders' = cpure_POP (Proxy::Proxy (Find ts)) (Fn $ Comp . (I <$>) . runDecoder (index find decoders))
  f repr = do
    childReprs <- decodeParent params repr
    children <- sequence_SOP' (ap_SOP decoders' childReprs)
    return $ to children

makeDecoders :: forall p g a s. (Real a, Floating a, Floating (Vector.Vector a), Usable a, All (And (KnownSizes s) (KnownCode s)) g) =>
  Complete p g -> NP (DecodeParams a s) g -> Rec (Decoder a s) p -> Rec (Decoder a s) (p :++ g)
makeDecoders complete params prim = fix f where
  f decoders = rAppend prim $ np2rec (cliftA2_NP (Proxy::Proxy (And (KnownSizes s) (KnownCode s))) (decodeRec decoders) (unAll_NP complete) params)

newtype AnyDecoder a s ts = AnyDecoder { runAnyDecoder :: forall t. Find ts t => ReprT a s t -> IO t }

makeDecoder :: (Real a, Floating a, Floating (Vector.Vector a), Usable a, All (And (KnownSizes s) (KnownCode s)) g) =>
  Complete p g -> NP (DecodeParams a s) g -> Rec (Decoder a s) p -> AnyDecoder a s (p :++ g)
makeDecoder complete params prim = AnyDecoder $ runDecoder (index find decoders)
  where decoders = makeDecoders complete params prim

newtype AutoDecoder a s t = AutoDecoder { runAutoDecoder :: Encoding a s t -> IO (Node a) }

autoDecodeRec :: forall ts a s t. (Floating a, Floating (Vector.Vector a), Usable a, KnownCode s t, All2 (KnownSize s) (Code t)) =>
  Rec (AutoDecoder a s) ts -> Dict (Contained ts) t -> DecodeParams a s t -> AutoDecoder a s t

autoDecodeRec autoDecoders Dict (DecodeParams aff params) = AutoDecoder f where
  autoDecoders' = cpure_POP (Proxy::Proxy (Find ts)) (Fn $ Comp . (K <$>) . runAutoDecoder (index find autoDecoders))
  f :: Encoding a s t -> IO (Node a)
  f (Generic parent children) = do
    Repr node <- affine' aff (Repr $ runReprT parent)
    let i = ns2int (unSOP children)

    prob <- makeSelect i node
    log_score <- makeUnary (negate . log) prob

    let childReprs = liftA_SOP getRepr children
    let decodings = cliftA_POP (Proxy::Proxy (KnownSize s)) (Comp . flip affineIn parent) params

    let dist :: forall t'. KnownSize s t' => (IO :.: ReprT a s) t' -> ReprT a s t' -> (IO :.: K (Node a)) t'
        dist dIO (ReprT c) = Comp $ do
          ReprT d <- unComp dIO
          diff <- makeBinary (-) d c
          norm <- makeDot diff diff
          return $ K norm

    norms <- collapse_SOP <$> (sequence_SOP' $ cliftA2_SOP (Proxy::Proxy (KnownSize s)) dist decodings childReprs)

    childNodes <- collapse_SOP <$> sequence_SOP' (ap_SOP autoDecoders' children)

    foldM (makeBinary (+)) log_score (norms ++ childNodes)

primAutoDecoder :: Num a => AutoDecoder a s t
primAutoDecoder = AutoDecoder f where
  f (Primitive _) = makeSource 0

makeAutoDecoders :: forall proxy p g a s. (Real a, Floating a, Floating (Vector.Vector a), Usable a, All (And (KnownSizes s) (KnownCode s)) g) =>
  Complete p g -> NP (DecodeParams a s) g -> Rec (AutoDecoder a s) p -> Rec (AutoDecoder a s) (p :++ g)
makeAutoDecoders complete params prim = fix f where
  --prim = np2rec $ pure_NP primitiveAutoDecoder :: Rec _ p
  f autoDecoders = rAppend prim $ np2rec (cliftA2_NP (Proxy::Proxy (And (KnownSizes s) (KnownCode s))) (autoDecodeRec autoDecoders) (unAll_NP complete) params)

newtype AnyAutoDecoder a s ts = AnyAutoDecoder { runAnyAutoDecoder :: forall t. Find ts t => Encoding a s t -> IO (Node a) }

makeAutoDecoder :: (Real a, Floating a, Floating (Vector.Vector a), Usable a, All (And (KnownSizes s) (KnownCode s)) g) =>
  Complete p g -> NP (DecodeParams a s) g -> Rec (AutoDecoder a s) p -> AnyAutoDecoder a s (p :++ g)
makeAutoDecoder complete params prim = AnyAutoDecoder $ runAutoDecoder (index find autoDecoders)
  where autoDecoders = makeAutoDecoders complete params prim

newtype AutoEncoder a s t = AutoEncoder { runAutoEncoder :: t -> IO (Node a) }

makeAutoEncoders :: forall p g a s. (Real a, Floating a, Floating (Vector.Vector a), Usable a, All (KnownSize s) g, All (And (KnownSizes s) (KnownCode s)) g) =>
  Complete p g -> NP (EncodeParams a s) g -> Rec (Encoder a s) p -> NP (DecodeParams a s) g -> Rec (AutoDecoder a s) p -> Rec (AutoEncoder a s) (p :++ g)

makeAutoEncoders complete encodeParams primEncoders decodeParams primAutoDecoders = autoEncoders where
  encoders = makeEncoders complete encodeParams primEncoders
  autoDecoders = makeAutoDecoders complete decodeParams primAutoDecoders

  makeAutoEncoder (Encoder e) (AutoDecoder d) = AutoEncoder f where
    f t = e t >>= d

  autoEncoders = rZipWith makeAutoEncoder encoders autoDecoders

newtype AnyAutoEncoder a s ts = AnyAutoEncoder { runAnyAutoEncoder :: forall t. Find ts t => t -> IO (Node a) }

makeAutoEncoder :: (Real a, Floating a, Floating (Vector.Vector a), Usable a, All (KnownSize s) g, All (And (KnownSizes s) (KnownCode s)) g) =>
  Complete p g -> NP (EncodeParams a s) g -> Rec (Encoder a s) p -> NP (DecodeParams a s) g -> Rec (AutoDecoder a s) p -> AnyAutoEncoder a s (p :++ g)
makeAutoEncoder complete encodeParams primEncoders decodeParams primAutoDecoders = AnyAutoEncoder $ runAutoEncoder (index find autoEncoders)
  where autoEncoders = makeAutoEncoders complete encodeParams primEncoders decodeParams primAutoDecoders


text:
{-# LANGUAGE ForeignFunctionInterface #-}
{-# CFILES cbits/vmath-zomplex.c #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Foreign.VMath.Zomplex
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--

module Foreign.VMath.Zomplex
    where

import Data.Complex ( Complex )
import Foreign.Ptr  ( Ptr )

foreign import ccall unsafe "vzConj"
    vzConj :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzShift"
    vzShift :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzAdd"
    vzAdd :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzSub"
    vzSub :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzAxpby"
    vzAxpby :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzScale"
    vzScale :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzMul"
    vzMul :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzNeg"
    vzNeg :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzAbs"
    vzAbs :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzSgn"
    vzSgn :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzInv"
    vzInv :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzDiv"
    vzDiv :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzExp"
    vzExp :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzSqrt"
    vzSqrt :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzLog"
    vzLog :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzPow"
    vzPow :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzSin"
    vzSin :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzCos"
    vzCos :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzTan"
    vzTan :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzASin"
    vzASin :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzACos"
    vzACos :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzATan"
    vzATan :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzSinh"
    vzSinh :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzCosh"
    vzCosh :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzTanh"
    vzTanh :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzASinh"
    vzASinh :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzACosh"
    vzACosh :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

foreign import ccall unsafe "vzATanh"
    vzATanh :: Int -> Ptr (Complex Double) -> Ptr (Complex Double) -> IO ()

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-|
Module      : MachineLearning.Utils.Report
Description : TIR expression data structures
Copyright   : (c) Fabricio Olivetti de Franca, 2022
License     : GPL-3
Maintainer  : fabricio.olivetti@gmail.com
Stability   : experimental
Portability : POSIX

Report generation.
-}
module MachineLearning.Utils.Report where

import           System.Directory
import           System.IO
import           System.Clock
import           Data.Maybe                  (fromMaybe)
import           Data.List                   (intercalate, foldl')
import           Data.List.Split             (splitOn)
import           Data.Int                    (Int64)
import qualified Data.Vector           as V
import qualified Data.Vector.Storable  as VS
import qualified Numeric.LinearAlgebra as LA

import Control.Evolution (Population)
import MachineLearning.Model.Measure                         (Measure(..))
import MachineLearning.TIR                         (Individual(..), assembleTree, replaceConsts)
import Data.SRTree.Print                            (showDefault, showPython)
import MachineLearning.Utils.Config                          (Config(..), Output(..), IOCfg(..), Task(..), AlgorithmCfg(..), getLogType, getMeasures)
import Data.SRTree (SRTree(..))

-- | Creates a file if it does not exist
createIfDoesNotExist :: FilePath -> IO Handle
createIfDoesNotExist fname = do
  isCreated <- doesFileExist fname
  if   isCreated
  then openFile fname AppendMode
  else openFile fname WriteMode

-- | writes all the stats from the final champion solution
writeChampionStats :: Config -> (Individual -> Maybe [Double]) -> Int64 -> Individual -> IO ()
writeChampionStats cfg fitTest totTime champion = do
  let dirname = case getLogType cfg of
                  Screen         -> "results/"
                  PartialLog dir -> dir
                  EvoLog dir     -> dir
      statsFname = dirname ++ "/stats.json"
      exprFname  = dirname ++ "/exprs.csv"
      cfgFname   = dirname ++ "/parameters.csv"
      task       = _task $ _algorithmCfg cfg
  createDirectoryIfMissing True dirname
  writeCfg cfgFname cfg
  writeExprs task exprFname champion
  writeStats statsFname fitTest (getMeasures cfg) totTime champion

-- | writes the stats to a file
writeStats :: FilePath -> (Individual -> Maybe [Double]) -> [Measure] -> Int64 -> Individual -> IO ()
writeStats statsFname fitTest measures totTime champion = do
    h <- createIfDoesNotExist statsFname
    hPutStrLn h $ "{" ++ intercalate "," json ++ "}"
    hClose h
  where
    mNames     = map _name measures
    trainNames = map (++"_train") mNames
    trainErrors = map show $ _fit champion
    testErrors  = map show $ fromMaybe nans $ fitTest champion
    testNames  = map (++"_test")  mNames
    fields     = "tot_time" : (trainNames ++ testNames)
    nans       = replicate (length testNames) (1/0)
    errors     = show totTime : (trainErrors ++ testErrors)
    toJson k v = show k ++ " : " ++ v
    json       = zipWith toJson fields errors

-- | creates a report file
writeCfg :: FilePath -> Config -> IO ()
writeCfg cfgFname cfg = do
  h <- createIfDoesNotExist cfgFname
  hPrint h cfg
  hClose h

-- | pretty write the expressions
writeExprs :: Task -> FilePath -> Individual -> IO ()
writeExprs task exprFname champion = do
    h <- createIfDoesNotExist exprFname
    hPutStrLn h $ unlines logs
    hPutStrLn h "======================================"
    hClose h
  where
    tir    = _chromo champion
    trees  = map getTree ws
    ws     = _weights champion
    bias   = V.head $ VS.convert $ head ws
    logs   = zipWith (\t w -> "\"" ++ showDefault t ++ "\",\"" ++ show w ++ "\",\"" ++ showPython t ++ "\"") trees ws

    getTree :: LA.Vector Double -> SRTree Int Double
    getTree w = let bias   = V.head $ VS.convert w
                    consts = V.tail $ VS.convert w
                    sigm z = 1 / (1+exp(-z))
                in  case task of
                      Classification _ -> sigm $ assembleTree bias $ replaceConsts tir consts
                      ClassMult      _ -> sigm $ assembleTree bias $ replaceConsts tir consts
                      _                -> assembleTree bias $ replaceConsts tir consts

-- | creates a log of the evolution process
evoLog :: Handle -> (Individual -> Maybe [Double]) -> Population Individual -> IO ()
evoLog h fitness pop = do
  let fitTrain   = V.toList $ V.map (head._fit) pop
      fitTest    = V.toList $ V.map (replaceWithNan . fitness) pop
      statsTrain = getStats fitTrain
      statsTest  = getStats fitTest
  hPutStrLn h $ statsToStr statsTrain ++ "," ++ statsToStr statsTest

getStats :: [Double] -> [Double]
getStats = postAgg . foldl' aggregate []

statsToStr :: Show a => [a] -> String
statsToStr = intercalate "," . map show

replaceWithNan :: Maybe [Double] -> Double
replaceWithNan Nothing      = 1/0
replaceWithNan (Just [])    = 1/0
replaceWithNan (Just (x:_)) = x

-- | Opens the first file available in the format "name.{i}.csv"
-- where 'i' follows a sequence from 0 onward.
openNext :: String -> IO Handle
openNext fname = go [fname ++ "." ++ show n ++ ".csv" | n <- [0 :: Int ..]]
  where
    -- this is a partial function applied to an infinite list
    -- so, what harm can it do?
    go []       = error "end of inifinity stream"
    go (fn:fns) = do b <- doesFileExist fn
                     if b
                        then go fns
                        else openFile fn WriteMode

postAgg :: [Double] -> [Double]
postAgg [best, worst, tot, count] = [best, worst, tot/count]
postAgg _ = error "wrong parameters count"

-- | aggregates the best, worst and average solutions
aggregate :: [Double] -> Double -> [Double]
aggregate [] train = [train,train,train,1]
aggregate [best, worst, tot, count] train = [min best train, max worst train, tot+train, count+1]
aggregate _ _ = error "wrong parameters count in aggregate"

-- | creates a log name instead of overwritting last one
createLogName :: Config -> FilePath
createLogName cfg = fileDir ++ "/" ++ dataset ++ "_evo"
  where
    fileDir = case getLogType cfg of
                EvoLog dir -> dir
                _           -> ""
    dataset = last $ splitOn "/" $ (_trainFilename . _ioCfg) cfg

-- | creates an IO logger function
makeLogger :: Config -> (Individual -> Maybe [Double]) -> IO (Population Individual -> IO (), Maybe Handle)
makeLogger cfg fitTest = case getLogType cfg of
                           EvoLog dir -> do createDirectoryIfMissing True dir
                                            h <- openNext (createLogName cfg)
                                            hPutStrLn h "bestTrain,worstTrain,avgTrain,bestTest,worstTest,avgTest"
                                            return (evoLog h fitTest, Just h)
                           _          -> return (\_ -> return (), Nothing)

-- | closes a file
closeIfJust :: Maybe Handle -> IO ()
closeIfJust Nothing  = return ()
closeIfJust (Just h) = hClose h

max_stars_repo_name: folivetti/tir

max_issues_repo_name: folivetti/tir

max_forks_repo_name: folivetti/tir


text:
{-|

 Module      :  HF
 Copyright   :  Copyright (c) David Schlegel
 License     :  BSD
 Maintainer  :  David Schlegel
 Stability   :  experimental
 Portability :  Haskell

 Module "HF".
-}

module HF (
{-| Module "HF.Data" provides all important datastructures and corresponding functions as well as dictionary lists.	-}
module HF.Data,
{-|In "ReadWrite" important functions are provided to

	* Read input (Geometry information, Basis set information) and convert to internal datastructures, provided in "Data"
	* Write output (Geometry, structure and basis set information)
-}
module HF.ReadWrite,
{-|Gaussian Integral evaluation plays an important role in quantum chemistry. In "HF.Gauss", functions will be provided to compute the most important integrals involving gaussian-type orbitals (GTOs). -}
module HF.Gauss)where
{-
--Total Packages included
import System.Process
import System.IO
import System.IO.Unsafe
import Data.List
import Data.List.Split
import Numeric.Container hiding (linspace)
import Data.Maybe
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data hiding (linspace)
import Text.Printf
-}

--Some Dictionaries
--dict1 = [("CARBON", "C"), ("NITROGEN", "N"), ("OXYGEN", "O"), ("HYDROGEN", "H")]
--dict2 = [("CARBON", 6.0), ("NITROGEN", 7.0), ("OXYGEN", 8.0), ("HYDROGEN", 1.0)]
--dict3 = [("S", 0), ("P", 1), ("D", 2), ("F", 3)]

import HF.Data

import HF.ReadWrite

import HF.Gauss

--import Matrices
import Numeric.Container hiding (linspace)

max_stars_repo_name: davidschlegel/hartree-fock

max_issues_repo_name: davidschlegel/hartree-fock

max_forks_repo_name: davidschlegel/hartree-fock


text:
{-# LANGUAGE TypeFamilies, BangPatterns, TypeOperators, FlexibleContexts, FlexibleInstances, ScopedTypeVariables #-}
module Main where

import Data.NeuralNetwork
import Data.NeuralNetwork.Backend.HMatrix
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Devel
import Parser

import Data.List(foldl',partition,maximumBy)
import Text.Printf (printf)
import Control.Monad
import Control.Monad.Except
import Text.PrettyPrint.Free hiding (flatten, (<>))
import System.Environment

import Data.Constraint
foo :: Backend b s => b -> s -> Inp (ConvertFromSpec b s) -> Env b (Out (ConvertFromSpec b s))
foo b s i = case witness b s of Dict -> compile b s >>= run . flip forward i

main = do x <- runExceptT $ do
                cnn <- compile ByHmatrix (In2D 28 28 :++ Reshape2DAs1D :++ FullConnect 256 :++ FullConnect 32 :++ FullConnect 10)
                debug cnn -- dotrain cnn >>= dotest
          case x of
            Left  _ -> putStrLn "Error."
            Right _ -> putStrLn "Done."

debug :: (Component n, Inp n ~ Image, Out n ~ Label, MonadIO m, RunInEnv (Run n) m)
      => n -> m ()
debug nn = do
    a0:a1:_ <- liftIO $ getArgs
    let cycle = read a0 :: Int
        rate  = read a1 :: Float
    liftIO $ putStrLn "Load training data."
    dataset <- liftIO $ uncurry zip <$> trainingData
    nn <- iterateM cycle (online rate dataset) nn
    liftIO $ putStrLn "Load test data."
    testset <- liftIO $ (uncurry zip <$> testData)
    flip mapM_ (take 10 testset) $ \(ds,ev) -> do
      pv <- run $ forward nn ds
      liftIO $ putStrLn ("+" ++ showPretty (text (printf "%02d:" $ postprocess pv) <+> pretty pv))
      liftIO $ putStrLn ("*" ++ showPretty (text (printf "%02d:" $ postprocess ev) <+> pretty ev))

dotrain :: (Component n, Inp n ~ Image, Out n ~ Label, MonadIO m, RunInEnv (Run n) m)
       => n -> m n
dotrain nn = do
    liftIO $ putStrLn "Load training data."
    dataset <- liftIO $ uncurry zip <$> trainingData
    liftIO $ putStrLn "Learning."
    iterateM 5 (online 0.0001 dataset) nn

dotest :: (Component n, Inp n ~ Image, Out n ~ Label, MonadIO m, RunInEnv (Run n) m)
       => n -> m ()
dotest !nn = do
    liftIO $ putStrLn "Load test data."
    testset <- liftIO $ (uncurry zip <$> testData)
    liftIO $ putStrLn "Start test"
    result  <- mapM ((postprocess <$>) . run . forward nn . fst) testset
    let expect = map (postprocess . snd) testset
        (co,wr)= partition (uncurry (==)) $ zip result expect
    liftIO $ putStrLn $ printf "correct: %d, wrong: %d" (length co) (length wr)

online :: (Component n, Inp n ~ Image, Out n ~ Label, MonadIO m, RunInEnv (Run n) m)
       => Float -> [(Inp n, Out n)] -> n -> m n
online rate ds !nn = walk ds nn
  where
    walk []     !nn = return nn
    walk (d:ds) !nn = do !nn <- run $ learn outcost' rate nn d
                         walk ds nn
    outcost' a b = return $ zipVectorWith cost' a b

iterateM :: (MonadIO m) => Int -> (a -> m a) -> a -> m a
iterateM n f x = walk 0 x
  where
    walk !i !a | i == n    = return a
               | otherwise = do -- when (i `mod` 10 == 0) $ putStrLn ("Iteration " ++ show i)
                                liftIO $ putStrLn ("Iteration " ++ show i)
                                a <- f a
                                walk (i+1) a

postprocess :: Vector Float -> Int
postprocess = fst . maximumBy cmp . zip [0..] . toList
  where cmp a b = compare (snd a) (snd b)

instance Pretty (Vector Float) where
  pretty vec = encloseSep langle rangle comma $ map (text . printf "%.04f") $ toList vec

showPretty x = displayS (renderPretty 0.4 500 x) ""

max_stars_repo_name: pierric/neural-network

max_issues_repo_name: pierric/neural-network

max_forks_repo_name: pierric/neural-network


text:
{-# OPTIONS_GHC -fno-warn-orphans #-}
{-# LANGUAGE CPP #-}
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
{-# LANGUAGE Trustworthy #-}
#endif
-----------------------------------------------------------------------------
-- |
-- Copyright   :  (C) 2012 Edward Kmett
-- License     :  BSD-style (see the file LICENSE)
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  provisional
-- Portability :  portable
--
-- Orphans
-----------------------------------------------------------------------------
module Linear.Instances () where

import Control.Applicative
import Control.Monad.Fix
import Control.Monad.Zip
import Data.Complex
import Data.Foldable
import Data.Functor.Bind
import Data.HashMap.Lazy as HashMap
import Data.Hashable
import Data.Semigroup
import Data.Semigroup.Foldable
import Data.Semigroup.Traversable
import Data.Traversable

instance (Hashable k, Eq k) => Apply (HashMap k) where
  (<.>) = HashMap.intersectionWith id

instance (Hashable k, Eq k) => Bind (HashMap k) where
  -- this is needlessly painful
  m >>- f = HashMap.fromList $ do
    (k, a) <- HashMap.toList m
    case HashMap.lookup k (f a) of
      Just b -> [(k,b)]
      Nothing -> []

instance Functor Complex where
  fmap f (a :+ b) = f a :+ f b
  {-# INLINE fmap #-}

instance Apply Complex where
  (a :+ b) <.> (c :+ d) = a c :+ b d

instance Applicative Complex where
  pure a = a :+ a
  (a :+ b) <*> (c :+ d) = a c :+ b d

instance Bind Complex where
  (a :+ b) >>- f = a' :+ b' where
    a' :+ _  = f a
    _  :+ b' = f b
  {-# INLINE (>>-) #-}

instance Monad Complex where
  return a = a :+ a
  {-# INLINE return #-}

  (a :+ b) >>= f = a' :+ b' where
    a' :+ _  = f a
    _  :+ b' = f b
  {-# INLINE (>>=) #-}

instance MonadZip Complex where
  mzipWith = liftA2

instance MonadFix Complex where
  mfix f = (let a :+ _ = f a in a) :+ (let _ :+ a = f a in a)

instance Foldable Complex where
  foldMap f (a :+ b) = f a `mappend` f b
  {-# INLINE foldMap #-}

instance Traversable Complex where
  traverse f (a :+ b) = (:+) <$> f a <*> f b
  {-# INLINE traverse #-}

instance Foldable1 Complex where
  foldMap1 f (a :+ b) = f a <> f b
  {-# INLINE foldMap1 #-}

instance Traversable1 Complex where
  traverse1 f (a :+ b) = (:+) <$> f a <.> f b
  {-# INLINE traverse1 #-}

max_stars_repo_name: abbradar/linear

max_issues_repo_name: abbradar/linear

max_forks_repo_name: abbradar/linear


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Graphics.Rendering.Plot
-- Copyright   :  (c) A. V. H. McPhail 2010
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- Graphical plots
--
-----------------------------------------------------------------------------

module Graphics.Rendering.Plot (
                                -- * Re-exported for convenience
                                module Graphics.Rendering.Plot.Figure.Simple
                               , module Graphics.Rendering.Plot.Figure
                               , module Graphics.Rendering.Plot.Render
                                -- * Example
                                -- $example
                               ) where

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

import Graphics.Rendering.Plot.Figure.Simple
import Graphics.Rendering.Plot.Figure
import Graphics.Rendering.Plot.Render

-----------------------------------------------------------------------------
{- $example

Create some data:

> ln = 25
> ts = linspace ln (0,1)
> rs = randomVector 0 Gaussian ln
>
> ss = sin (15*2*pi*ts)
> ds = 0.25*rs + ss
> es = constant (0.25*(stddev rs)) ln
>
> fs :: Double -> Double
> fs = sin . (15*2*pi*)

Perform actions in 'Figure a' to create a figure

> test_graph = do
>         withTextDefaults $ setFontFamily "OpenSymbol"
>         withTitle $ setText "Testing plot package:"
>         withSubTitle $ do
>                        setText "with 1 second of a 15Hz sine wave"
>                        setFontSize 10
>         setPlots 1 1
>         withPlot (1,1) $ do
>                          setDataset (ts,[point (ds,es) (Cross,red),line fs blue])
>                          addAxis XAxis (Side Lower) $ withAxisLabel $ setText "time (s)"
>                          addAxis YAxis (Side Lower) $ withAxisLabel $ setText "amplitude"
>                          addAxis XAxis (Value 0) $ return ()
>                          setRangeFromData XAxis Lower Linear
>                          setRange YAxis Lower Linear (-1.25) 1.25

Render the graph to a Cairo 'Render ()' action that takes the width
and height of the drawing area

> test_render :: (Double,Double) -> Render ()
> test_render = render test_graph

The same graph using the 'Simple' interface

> test_graph2 = do
>         plot (ts,[point (ds,es) (Cross,red),line fs blue])
>         title "Testing plot package:"
>         subtitle "with 1 second of a 15Hz sine wave"
>         xlabel "time (s)"
>         ylabel "amplitude"
>         yrange Linear (-1.25) 1.25

The 'Render a' action can be used in GTK or with Cairo to write to file in PS, PDF, SVG, or PNG

Display a greyscale matrix

> ms :: Matrix Double
> ms = buildMatrix 64 64 (\(x,y) -> sin (2*2*pi*(fromIntegral x)/64) * cos (5*2*pi*(fromIntegral y)/64))

> mat_fig = do
>         setPlots 1 1
>         withPlot (1,1) $ do
>                          setDataset ms
>                          addAxis XAxis (Side Lower) $ setTickLabelFormat "%.0f"
>                          addAxis YAxis (Side Lower) $ setTickLabelFormat "%.0f"
>                          setRangeFromData XAxis Lower Linear
>                          setRangeFromData YAxis Lower Linear

The ODE example from hmatrix:

> import Numeric.GSL
> import Numeric.LinearAlgebra

> xdot t [x,v] = [v, -0.95*x - 0.1*v]
> ts = linspace 100 (0,20)
> sol = odeSolve xdot [10,0] ts

> ode_fig = plot (Line,ts,[sol])

-}

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

max_stars_repo_name: JackTheEngineer/plot


text:
module Statistics.Information.Discrete.MutualInfo where

import Data.Matrix
import Statistics.Information.Discrete.Entropy

mi :: Eq a => Int -> Matrix a -> Matrix a -> Double
mi base xs ys = entropy base xs + entropy base ys - entropy base (xs <|> ys)

cmi :: Eq a => Int -> Matrix a -> Matrix a -> Matrix a -> Double
cmi base xs ys zs = hxz + hyz - hxyz - hz where
  hxz = entropy base (xs <|> zs)
  hyz = entropy base (ys <|> zs)
  hxyz = entropy base (xs <|> ys <|> zs)
  hz = entropy base zs

max_stars_repo_name: eligottlieb/Shannon

max_issues_repo_name: eligottlieb/Shannon

max_forks_repo_name: eligottlieb/Shannon


text:
{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE GADTs                    #-}
module Grenade.Layers.Internal.Update (
    descendMatrix
  , descendVector
  , MatrixInputValues (..)
  , MatrixResult (..)
  , VectorInputValues (..)
  , VectorResult (..)
  ) where

import           Data.Maybe                   (fromJust)
import qualified Data.Vector.Storable         as U (unsafeFromForeignPtr0,
                                                    unsafeToForeignPtr0)
import           Foreign                      (mallocForeignPtrArray, withForeignPtr)
import           Foreign.Ptr                  (Ptr)
import           GHC.TypeLits
import           Numeric.LinearAlgebra        (Vector, flatten)
import qualified Numeric.LinearAlgebra.Devel  as U
import           Numeric.LinearAlgebra.Static
import           System.IO.Unsafe             (unsafePerformIO)

import           Grenade.Core.Optimizer
import           Grenade.Types

data MatrixInputValues rows columns
  = MatrixValuesSGD
      !(L rows columns) -- ^ current weights
      !(L rows columns) -- ^ gradients
      !(L rows columns) -- ^ last update (old momentum)
  | MatrixValuesAdam
      !Int   -- ^ Step
      !(L rows columns) -- ^ current weights
      !(L rows columns) -- ^ gradients
      !(L rows columns) -- ^ current m
      !(L rows columns) -- ^ current v

data MatrixResult rows columns
  = MatrixResultSGD
      { matrixActivations :: !(L rows columns) -- ^ new activations (weights)
      , matrixMomentum    :: !(L rows columns) -- ^ new momentum
      }
  | MatrixResultAdam
      { matrixActivations :: !(L rows columns) -- ^ new activations (weights)
      , matrixM           :: !(L rows columns) -- ^ new m
      , matrixV           :: !(L rows columns) -- ^ new v
      }

data VectorInputValues r
  = VectorValuesSGD
      !(R r) -- ^ current weights
      !(R r) -- ^ gradients
      !(R r) -- ^ last update (old momentum)
  | VectorValuesAdam
      !Int   -- ^ Step
      !(R r) -- ^ current weights
      !(R r) -- ^ current gradients
      !(R r) -- ^ current m
      !(R r) -- ^ current v

data VectorResult r
  = VectorResultSGD
      { vectorBias     :: !(R r) -- ^ new activations (bias)
      , vectorMomentum :: !(R r) -- ^ new momentum
      }
  | VectorResultAdam
      { vectorBias :: !(R r) -- ^ new activations (bias)
      , vectorM    :: !(R r) -- ^ new m
      , vectorV    :: !(R r) -- ^ new v
      }

descendMatrix :: (KnownNat rows, KnownNat columns) => Optimizer o -> MatrixInputValues rows columns -> MatrixResult rows columns
descendMatrix (OptSGD rate momentum regulariser) (MatrixValuesSGD weights gradient lastUpdate) =
  let (rows, cols) = size weights
      len          = rows * cols
      -- Most gradients come in in ColumnMajor,
      -- so we'll transpose here before flattening them
      -- into a vector to prevent a copy.
      --
      -- This gives ~15% speed improvement for LSTMs.
      weights'     = flatten . tr . extract $ weights
      gradient'    = flatten . tr . extract $ gradient
      lastUpdate'  = flatten . tr . extract $ lastUpdate
      (vw, vm)     = descendUnsafeSGD len rate momentum regulariser weights' gradient' lastUpdate'

      -- Note that it's ColumnMajor, as we did a transpose before
      -- using the internal vectors.
      mw           = U.matrixFromVector U.ColumnMajor rows cols vw
      mm           = U.matrixFromVector U.ColumnMajor rows cols vm
  in  MatrixResultSGD (fromJust . create $ mw) (fromJust . create $ mm)
descendMatrix (OptAdam alpha beta1 beta2 epsilon lambda) (MatrixValuesAdam step weights gradient m v) =
  let (rows, cols) = size weights
      len          = rows * cols
      -- Most gradients come in in ColumnMajor,
      -- so we'll transpose here before flattening them
      -- into a vector to prevent a copy.
      --
      -- This gives ~15% speed improvement for LSTMs.
      weights'  = flatten . tr . extract $ weights
      gradient' = flatten . tr . extract $ gradient
      m'        = flatten . tr . extract $ m
      v'        = flatten . tr . extract $ v
      (vw, vm, vv)     = descendUnsafeAdam len step alpha beta1 beta2 epsilon lambda weights' gradient' m' v'

      -- Note that it's ColumnMajor, as we did a transpose before
      -- using the internal vectors.
      mw           = U.matrixFromVector U.ColumnMajor rows cols vw
      mm           = U.matrixFromVector U.ColumnMajor rows cols vm
      mv           = U.matrixFromVector U.ColumnMajor rows cols vv
  in  MatrixResultAdam (fromJust . create $ mw) (fromJust . create $ mm) (fromJust . create $ mv)
descendMatrix opt _ = error $ "optimzer does not match to MatrixInputValues in implementation! Optimizer: " ++ show opt

descendVector :: (KnownNat r) => Optimizer o -> VectorInputValues r -> VectorResult r
descendVector (OptSGD rate momentum regulariser) (VectorValuesSGD weights gradient lastUpdate) =
  let len          = size weights
      weights'     = extract weights
      gradient'    = extract gradient
      lastUpdate'  = extract lastUpdate
      (vw, vm)     = descendUnsafeSGD len rate momentum regulariser weights' gradient' lastUpdate'
  in  VectorResultSGD (fromJust $ create vw) (fromJust $ create vm)
descendVector (OptAdam alpha beta1 beta2 epsilon lambda) (VectorValuesAdam step weights gradient m v) =
  let len       = size weights
      weights'  = extract weights
      gradient' = extract gradient
      m'        = extract m
      v'        = extract v
      (vw, vm, vv)     = descendUnsafeAdam len step alpha beta1 beta2 epsilon lambda weights' gradient' m' v'
  in  VectorResultAdam (fromJust $ create vw) (fromJust $ create vm) (fromJust $ create vv)
descendVector opt _ = error $ "optimzer does not match to VectorInputValues in implementation! Optimizer: " ++ show opt

descendUnsafeSGD :: Int -> RealNum -> RealNum -> RealNum -> Vector RealNum -> Vector RealNum -> Vector RealNum -> (Vector RealNum, Vector RealNum)
descendUnsafeSGD len rate momentum regulariser weights gradient lastUpdate =
  unsafePerformIO $ do
    outWPtr <- mallocForeignPtrArray len
    outMPtr <- mallocForeignPtrArray len
    let (wPtr, _) = U.unsafeToForeignPtr0 weights
    let (gPtr, _) = U.unsafeToForeignPtr0 gradient
    let (lPtr, _) = U.unsafeToForeignPtr0 lastUpdate

    withForeignPtr wPtr $ \wPtr' ->
      withForeignPtr gPtr $ \gPtr' ->
        withForeignPtr lPtr $ \lPtr' ->
          withForeignPtr outWPtr $ \outWPtr' ->
            withForeignPtr outMPtr $ \outMPtr' ->
              descend_sgd_cpu len rate momentum regulariser wPtr' gPtr' lPtr' outWPtr' outMPtr'

    return (U.unsafeFromForeignPtr0 outWPtr len, U.unsafeFromForeignPtr0 outMPtr len)

descendUnsafeAdam ::
     Int -- Len
  -> Int -- Step
  -> RealNum -- Alpha
  -> RealNum -- Beta1
  -> RealNum -- Beta2
  -> RealNum -- Epsilon
  -> RealNum -- Lambda
  -> Vector RealNum -- Weights
  -> Vector RealNum -- Gradient
  -> Vector RealNum -- M
  -> Vector RealNum -- V
  -> (Vector RealNum, Vector RealNum, Vector RealNum)
descendUnsafeAdam len step alpha beta1 beta2 epsilon lambda weights gradient m v =
  unsafePerformIO $ do
    outWPtr <- mallocForeignPtrArray len
    outMPtr <- mallocForeignPtrArray len
    outVPtr <- mallocForeignPtrArray len
    let (wPtr, _) = U.unsafeToForeignPtr0 weights
    let (gPtr, _) = U.unsafeToForeignPtr0 gradient
    let (mPtr, _) = U.unsafeToForeignPtr0 m
    let (vPtr, _) = U.unsafeToForeignPtr0 v
    withForeignPtr wPtr $ \wPtr' ->
      withForeignPtr gPtr $ \gPtr' ->
        withForeignPtr mPtr $ \mPtr' ->
          withForeignPtr vPtr $ \vPtr' ->
            withForeignPtr outWPtr $ \outWPtr' ->
              withForeignPtr outMPtr $ \outMPtr' ->
                withForeignPtr outVPtr $ \outVPtr' -> descend_adam_cpu len step alpha beta1 beta2 epsilon lambda wPtr' gPtr' mPtr' vPtr' outWPtr' outMPtr' outVPtr'
    return (U.unsafeFromForeignPtr0 outWPtr len, U.unsafeFromForeignPtr0 outMPtr len, U.unsafeFromForeignPtr0 outVPtr len)

foreign import ccall unsafe
    descend_sgd_cpu
      :: Int -> RealNum -> RealNum -> RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> IO ()

foreign import ccall unsafe
    descend_adam_cpu
      :: Int -> Int -> RealNum -> RealNum -> RealNum -> RealNum -> RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> IO ()

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE RankNTypes #-}
{-|
Module: MachineLearning.NeuralNetwork.Layer
Description: Neural Network's Layer
Copyright: (c) Alexander Ignatyev, 2017-2018.
License: BSD-3
Stability: experimental
Portability: POSIX

Neural Network's Layer
-}

module MachineLearning.NeuralNetwork.Layer
(
  Layer(..)
  , Cache(..)
  , affineForward
  , affineBackward
)

where

import Prelude hiding ((<>))
import MachineLearning.Types (R, Matrix)
import MachineLearning.Utils (sumByColumns)
import qualified Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra ((<>))
import qualified Control.Monad.Random as RndM

data Cache = Cache {
  cacheZ :: Matrix
  , cacheX :: Matrix
  , cacheW :: Matrix
  };

data Layer = Layer {
  lUnits :: Int
  , lForward :: Matrix -> Matrix -> Matrix -> Matrix
  , lBackward :: Matrix -> Cache -> (Matrix, Matrix, Matrix)
  , lActivation :: Matrix -> Matrix
  , lActivationGradient :: Matrix -> Matrix -> Matrix
  , lInitializeThetaM :: forall g. RndM.RandomGen g => (Int, Int) -> RndM.Rand g (Matrix, Matrix)
  }

affineForward :: Matrix -> Matrix -> Matrix -> Matrix
affineForward x b w = (x <> LA.tr w) + b

affineBackward delta (Cache _ x w) = (dx, db, dw)
  where m = fromIntegral $ LA.rows x
        dx = delta <> w
        db = (sumByColumns delta)/m
        dw = (LA.tr delta <> x)/m

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE BangPatterns, FlexibleContexts, UnboxedTuples #-}
-- |
-- Module    : Statistics.Sample.KernelDensity
-- Copyright : (c) 2011 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Kernel density estimation.  This module provides a fast, robust,
-- non-parametric way to estimate the probability density function of
-- a sample.
--
-- This estimator does not use the commonly employed \"Gaussian rule
-- of thumb\".  As a result, it outperforms many plug-in methods on
-- multimodal samples with widely separated modes.

module Statistics.Sample.KernelDensity
    (
    -- * Estimation functions
      kde
    -- , kde_
    -- * References
    -- $references
    ) where

import Numeric.MathFunctions.Constants (m_sqrt_2_pi)
import Prelude hiding (const, min, max, sum)
import Statistics.Function (minMax, nextHighestPowerOfTwo)
import Statistics.Math.RootFinding (fromRoot, ridders)
import Statistics.Sample.Histogram (histogram_)
import Statistics.Sample.Internal (sum)
import Statistics.Transform (CD, dct, idct)
import qualified Data.Vector.Generic  as G
import qualified Data.Vector.Unboxed  as U
import qualified Data.Vector          as V

-- | Gaussian kernel density estimator for one-dimensional data, using
-- the method of Botev et al.
--
-- The result is a pair of vectors, containing:
--
-- * The coordinates of each mesh point.  The mesh interval is chosen
--   to be 20% larger than the range of the sample.  (To specify the
--   mesh interval, use 'kde_'.)
--
-- * Density estimates at each mesh point.
kde :: (G.Vector v CD, G.Vector v Double, G.Vector v Int)
    => Int
    -- ^ The number of mesh points to use in the uniform discretization
    -- of the interval @(min,max)@.  If this value is not a power of
    -- two, then it is rounded up to the next power of two.
    -> v Double -> (v Double, v Double)
kde n0 xs = kde_ n0 (lo - range / 10) (hi + range / 10) xs
  where
    (lo,hi) = minMax xs
    range   | G.length xs <= 1 = 1       -- Unreasonable guess
            | lo == hi         = 1       -- All elements are equal
            | otherwise        = hi - lo
{-# INLINABLE  kde #-}
{-# SPECIAlIZE kde :: Int -> U.Vector Double -> (U.Vector Double, U.Vector Double) #-}
{-# SPECIAlIZE kde :: Int -> V.Vector Double -> (V.Vector Double, V.Vector Double) #-}

-- | Gaussian kernel density estimator for one-dimensional data, using
-- the method of Botev et al.
--
-- The result is a pair of vectors, containing:
--
-- * The coordinates of each mesh point.
--
-- * Density estimates at each mesh point.
kde_ :: (G.Vector v CD, G.Vector v Double, G.Vector v Int)
     => Int
     -- ^ The number of mesh points to use in the uniform discretization
     -- of the interval @(min,max)@.  If this value is not a power of
     -- two, then it is rounded up to the next power of two.
     -> Double
     -- ^ Lower bound (@min@) of the mesh range.
     -> Double
     -- ^ Upper bound (@max@) of the mesh range.
     -> v Double
     -> (v Double, v Double)
kde_ n0 min max xs
  | G.null xs = error "Statistics.KernelDensity.kde: empty sample"
  | n0 <= 1   = error "Statistics.KernelDensity.kde: invalid number of points"
  | otherwise = (mesh, density)
  where
    mesh = G.generate ni $ \z -> min + (d * fromIntegral z)
        where d = r / (n-1)
    density = G.map (/(2 * r)) . idct $ G.zipWith f a (G.enumFromTo 0 (n-1))
      where f b z = b * exp (sqr z * sqr pi * t_star * (-0.5))
    !n  = fromIntegral ni
    !ni = nextHighestPowerOfTwo n0
    !r  = max - min
    a   = dct . G.map (/ sum h) $ h
        where h = G.map (/ len) $ histogram_ ni min max xs
    !len    = fromIntegral (G.length xs)
    !t_star = fromRoot (0.28 * len ** (-0.4)) . ridders 1e-14 (0,0.1) $ \x ->
              x - (len * (2 * sqrt pi) * go 6 (f 7 x)) ** (-0.4)
      where
        f q t = 2 * pi ** (q*2) * sum (G.zipWith g iv a2v)
          where g i a2 = i ** q * a2 * exp ((-i) * sqr pi * t)
                a2v = G.map (sqr . (*0.5)) $ G.tail a
                iv = G.map sqr $ G.enumFromTo 1 (n-1)
        go s !h | s == 1    = h
                | otherwise = go (s-1) (f s time)
          where time  = (2 * const * k0 / len / h) ** (2 / (3 + 2 * s))
                const = (1 + 0.5 ** (s+0.5)) / 3
                k0    = U.product (G.enumFromThenTo 1 3 (2*s-1)) / m_sqrt_2_pi
    sqr x = x * x
{-# INLINABLE  kde_ #-}
{-# SPECIAlIZE kde_ :: Int -> Double -> Double -> U.Vector Double -> (U.Vector Double, U.Vector Double) #-}
{-# SPECIAlIZE kde_ :: Int -> Double -> Double -> V.Vector Double -> (V.Vector Double, V.Vector Double) #-}

-- $references
--
-- Botev. Z.I., Grotowski J.F., Kroese D.P. (2010). Kernel density
-- estimation via diffusion. /Annals of Statistics/
-- 38(5):2916–2957. <http://arxiv.org/pdf/1011.2602>

max_stars_repo_name: runeksvendsen/hs-gauge

max_issues_repo_name: runeksvendsen/hs-gauge

max_forks_repo_name: runeksvendsen/hs-gauge


text:
{-# LANGUAGE StandaloneDeriving, DeriveGeneric #-}

module T5884 where

import GHC.Generics

import Data.Complex

deriving instance Generic (Complex v)

max_stars_repo_name: ocharles/ghc

max_issues_repo_name: ocharles/ghc

max_forks_repo_name: ocharles/ghc


text:
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}

module HVX.BlackBoxTests.HuberTestSpec where

import Test.Hspec
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Util (zeros)

import HVX
import HVX.Internal.TestUtil

-- Problem definition.
n = 4
a = EConst $ (n><n)
    ([0.5377  , 0.3188  , 3.5784  , 0.7254
    , 1.8339  , -1.3077 , 2.7694  , -0.0631
    , -2.2588 , -0.4336 , -1.3499 , 0.7147
    , 0.8622  , 0.3426  , 3.0349  , -0.2050 ] :: [Double])
b = EConst $
    (n><1) ([ -1.0689 , -0.8095 , -2.9443 , -1.4384 ] :: [Double])
c = EConst $ (n><n)
    ([ -0.1241 , 0.6715  , 0.4889  , 0.2939
    , 1.4897  , -1.2075 , 1.0347  , -0.7873
    , 1.4090  , 0.7172  , 0.7269  , 0.8884
    , 1.4172  , 1.6302  , -0.3034 , -1.1471 ] :: [Double])
d = EConst $ (n><1)
    ([ 0.3252  , -0.7549 , 1.3703  , -1.7115 ] :: [Double])

x = EVar "x"
y = EVar "y"
constZeroVector = EConst $ zeros n 1

subgradAns = subgradMinimize
  (hmax $ huber 3 (a *~ x +~ b) +~ berhu 1 (c *~ y +~ d))
  [x +~ y <=~ constZeroVector]
  (decNonSumStep 1.0) 100000
  [("x", zeros n 1), ("y", zeros n 1)]

ellipsoidAns = ellipsoidMinimize
  (hmax $ huber 3 (a *~ x +~ b) +~ berhu 1 (c *~ y +~ d))
  [x +~ y <=~ constZeroVector]
  [("x", n), ("y", n)]
  1e-16 1e10

(subgradVars, subgradOptval) = subgradAns

(ellipsoidVars, ellipsoidOptval, ellipsoidUBound) = ellipsoidAns

-- CVX's results.
cvxOptval = 0.404882
--cvxOptx = (n><1) [-1.6137 ,-1.1202 ,0.7986 ,-0.7487 ]
--cvxOpty = (n><1) [0.2428 ,-0.0191 ,-0.7986 ,-1.0079 ]

-- Verify that HVX matches CVX.
spec :: Spec
spec =
  describe "Verify that HVX matches CVX for huber/berhu" $ do
    it "Verify that HVX subgrad matches CVX for huber/berhu" $
      subgradOptval `shouldSatisfy` fpequalsApprox cvxOptval
    it "Verify that HVX ellipsoid matches CVX for huber/berhu" $
      ellipsoidOptval `shouldSatisfy` fpequalsApprox cvxOptval

main :: IO ()
main = hspec spec

max_stars_repo_name: wellposed/hvx

max_issues_repo_name: wellposed/hvx

max_forks_repo_name: wellposed/hvx


text:
module Numeric.Physics.Particle where

import Numeric.LinearAlgebra.Vector

data Particle a = Particle
  { pPos     :: !(Vec3 a) -- ^ Position of particle in the world
  , pVel     :: !(Vec3 a) -- ^ Velocity of particle
  , pAcc     :: !(Vec3 a) -- ^ Acceleration of particle
  , pDamp    :: !a        -- ^ damping factor for this particle
  , pInvMass :: !a        -- ^ inverse mass (1/m) of this particle
  } deriving (Read,Show,Ord,Eq)

-- | Takes a particle, a force, and a time duration
-- and gives the particle updated for time slice
integrate :: Floating a => Particle a -> Vec3 a -> a -> Particle a
integrate p f dt =
  let pos  = (dt         *> pVel p) <+> pPos p
      acc' = (pInvMass p *> f)      <+> pAcc p
      vel  = (dt         *> acc')   <+> pVel p
  in p { pPos = pos
       , pVel = ((pDamp p)**dt) *> vel
       }
{-# SPECIALIZE INLINE integrate :: Particle Double -> Vec3 Double -> Double -> Particle Double #-}
{-# SPECIALIZE INLINE integrate :: Particle Float  -> Vec3 Float  -> Float  -> Particle Float  #-}

max_stars_repo_name: dagit/physics

max_issues_repo_name: dagit/physics

max_forks_repo_name: dagit/physics


text:
{-# LANGUAGE RecordWildCards #-}

module Trade.Algorithm.Drawdown where

import qualified Data.List as List
import qualified Data.Vector as Vec

import qualified Statistics.Sample as Sample

{-

dds :: [Double] -> [(Int, Double)]
dds xs =
  let go n a [] = [(n, a)]
      go n 0 (x:xs) = go n x xs
      go n a (0:xs) = go (n+1) a xs
      go n a (x:xs) =
        case a*x > 0 of
          True -> go (n+1) (a+x) xs
          False -> let ys = go 1 x xs in (n, a):ys
  in go 1 0 (zipWith (-) (tail xs) xs)

drawdownSeries :: [Double] -> [(Double, Double)]
drawdownSeries us@(u:_) =
  let go acc [] = [acc]
      go (a, b) ((x, y):us) =
        let vs = go (a+fromIntegral x, b+y) us
        in (a, b) : vs
  in go (0, u) (dds us)

drawdowns :: [Double] -> [(Int, Double)]
drawdowns = filter ((<0) . snd) . dds

numberOfDrawdowns :: [Double] -> Int
numberOfDrawdowns = length . drawdowns

maximumDrawdown :: [Double] -> Double
maximumDrawdown = minimum . map snd . dds

wins :: [Double] -> [(Int, Double)]
wins = filter ((>0) . snd) . dds

numberOfWins :: [Double] -> Int
numberOfWins = length . wins

maximumWin :: [Double] -> Double
maximumWin = maximum . map snd . dds

data DrawdownStatistics = DrawdownStatistics {
  numOfWins :: !Int
  , maxWin :: !Double
  , meanWin :: !Double
  , stdDevWin :: !Double
  , numOfDrawdowns :: !Int
  , maxDrawdown :: !Double
  , meanDrawdown :: !Double
  , stdDevDrawdown :: !Double
  } deriving (Show)

drawdownStats :: [Double] -> DrawdownStatistics
drawdownStats xs =
  let ds = map snd (dds xs)
      wins = filter (>0) ds
      vwins = Vec.fromList wins
      losses = filter (<0) ds
      vlosses = Vec.fromList losses
  in DrawdownStatistics {
    numOfWins = length wins
    , maxWin = maximum wins
    , meanWin = Sample.mean vwins
    , stdDevWin = Sample.stdDev vwins
    , numOfDrawdowns = length losses
    , maxDrawdown = minimum losses
    , meanDrawdown = Sample.mean vlosses
    , stdDevDrawdown = Sample.stdDev vlosses
    }

reportDrawdownStats :: [Double] -> HtmlIO
reportDrawdownStats xs =
  let DrawdownStatistics{..} = drawdownStats xs
  in Rep.vtable $
     ["Number of wins", show numOfWins]
     : ["Win max", show maxWin]
     : ["Win mean", show meanWin]
     : ["Win stddev", show stdDevWin]
     : ["----------"]
     : ["Number of drawdowns", show numOfDrawdowns]
     : ["Drawdown max", show maxDrawdown]
     : ["Drawdown mean", show meanDrawdown]
     : ["Drawdown stddev", show stdDevDrawdown]
     : []

-}

max_stars_repo_name: fphh/trade

max_issues_repo_name: fphh/trade

max_forks_repo_name: fphh/trade


text:
{-# LANGUAGE GADTs #-}

module HVX.Internal.SymbolicSubgrad where

import Numeric.LinearAlgebra hiding (i)
import Numeric.LinearAlgebra.Util

import HVX.Internal.Matrix
import HVX.Internal.Primitives

type Vars = [(Var, Mat)]

evaluate :: Expr vex mon -> Vars -> Mat
evaluate (EConst c) _ = c
evaluate (EVar var) vars
    | Just val <- lookup var vars = val
    | otherwise = error $ "Variable " ++ show var ++ " not set."
evaluate (EFun f e) vars = getFun f (evaluate e vars)
evaluate (EAdd l r) vars = evaluate l vars + evaluate r vars

jacobianWrtVar :: Expr vex mon  -- ^ Expression whose jacobian to take.
  -> Vars                       -- ^ Variables appearing in the expression.
  -> Var                        -- ^ Variable to take jacobian with respect to.
  -> Mat                        -- ^ The jacobian.
jacobianWrtVar (EConst c) vars wrtVar
    | Just wrtVal <- lookup wrtVar vars = zeros (rows c) (rows wrtVal)
    | otherwise = error $ "Variable " ++ show wrtVar ++ " not set."
jacobianWrtVar (EVar var) vars wrtVar
    | var == wrtVar, Just val <- lookup var vars = ident $ rows val
    | Just val <- lookup var vars, Just wrtVal <- lookup wrtVar vars = zeros (rows val) (rows wrtVal)
    | otherwise = error $ "Variable " ++ show var ++ " or " ++ show wrtVar ++ " not set."
-- Chain rule: ddx f(e) = dde f * ddx e
jacobianWrtVar (EFun f e) vars var = dde_f <> ddx_e
  where
    dde_f = getJacobian f val
    ddx_e = jacobianWrtVar e vars var
    val   = evaluate e vars
-- Special case for addition.
jacobianWrtVar (EAdd l r) vars var = ddx_l + ddx_r
  where
    ddx_l = jacobianWrtVar l vars var
    ddx_r = jacobianWrtVar r vars var