text:
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleContexts      #-}
{-|
Module      : Grenade.Layers.Reshape
Description : Multipurpose reshaping layer
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Reshape (
    Reshape (..)
  ) where

import           Data.Serialize

import           Data.Singletons.TypeLits
import           GHC.TypeLits
import           Control.DeepSeq              (NFData (..))
import           GHC.Generics                 (Generic)

import           Numeric.LinearAlgebra.Static
import           Numeric.LinearAlgebra.Data as LA ( flatten )

import           Grenade.Core

-- | Reshape Layer
--
-- The Reshape layer can flatten any 2D or 3D image to 1D vector with the
-- same number of activations, as well as cast up from 1D to a 2D or 3D
-- shape.
--
-- Can also be used to turn a 3D image with only one channel into a 2D image
-- or vice versa.
data Reshape = Reshape
  deriving (Show,Generic,NFData)

instance UpdateLayer Reshape where
  type Gradient Reshape = ()
  runUpdate _ _ _ = Reshape

instance RandomLayer Reshape where
  createRandomWith _ _ = return Reshape

instance (KnownNat a, KnownNat x, KnownNat y, a ~ (x * y)) => Layer Reshape ('D2 x y) ('D1 a) where
  type Tape Reshape ('D2 x y) ('D1 a) = ()
  runForwards _ (S2D y)   =  ((), fromJust' . fromStorable . flatten . extract $ y)
  runBackwards _ _ (S1D y) = ((), fromJust' . fromStorable . extract $ y)

instance (KnownNat a, KnownNat x, KnownNat y, KnownNat (x * z), KnownNat z, a ~ (x * y * z)) => Layer Reshape ('D3 x y z) ('D1 a) where
  type Tape Reshape ('D3 x y z) ('D1 a) = ()
  runForwards _ (S3D y)     = ((), fromJust' . fromStorable . flatten . extract $ y)
  runBackwards _ _ (S1D y) = ((), fromJust' . fromStorable . extract $ y)

instance (KnownNat y, KnownNat x, KnownNat z, z ~ 1) => Layer Reshape ('D3 x y z) ('D2 x y) where
  type Tape Reshape ('D3 x y z) ('D2 x y) = ()
  runForwards _ (S3D y)    = ((), S2D y)
  runBackwards _ _ (S2D y) = ((), S3D y)

instance (KnownNat y, KnownNat x, KnownNat z, z ~ 1) => Layer Reshape ('D2 x y) ('D3 x y z) where
  type Tape Reshape ('D2 x y) ('D3 x y z) = ()
  runForwards _ (S2D y)    = ((), S3D y)
  runBackwards _ _ (S3D y) = ((), S2D y)

instance (KnownNat a, KnownNat x, KnownNat y, a ~ (x * y)) => Layer Reshape ('D1 a) ('D2 x y) where
  type Tape Reshape ('D1 a) ('D2 x y) = ()
  runForwards _ (S1D y)   =  ((), fromJust' . fromStorable . extract $ y)
  runBackwards _ _ (S2D y) = ((), fromJust' . fromStorable . flatten . extract $ y)

instance (KnownNat a, KnownNat x, KnownNat y, KnownNat (x * z), KnownNat z, a ~ (x * y * z)) => Layer Reshape ('D1 a) ('D3 x y z) where
  type Tape Reshape ('D1 a) ('D3 x y z) = ()
  runForwards _ (S1D y)     = ((), fromJust' . fromStorable . extract $ y)
  runBackwards _ _ (S3D y)  = ((), fromJust' . fromStorable . flatten . extract $ y)

instance Serialize Reshape where
  put _ = return ()
  get = return Reshape

fromJust' :: Maybe x -> x
fromJust' (Just x) = x
fromJust' Nothing  = error $ "Reshape error: data shape couldn't be converted."

-------------------- GNum instances --------------------

instance GNum Reshape where
  _ |* Reshape = Reshape
  _ |+ Reshape = Reshape
  gFromRational _ = Reshape

max_stars_repo_name: koenigmaximilian/grenade

max_issues_repo_name: koenigmaximilian/grenade

max_forks_repo_name: koenigmaximilian/grenade


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Types
-- Copyright : (c) 2009 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Data types common used in statistics
module Statistics.Types
    ( -- * Confidence level
      CL
      -- ** Accessors
    , confidenceLevel
    , significanceLevel
      -- ** Constructors
    , mkCL
      -- ** Constants and conversion to nσ
    , cl95
      -- * Estimates and upper/lower limits
    , Estimate(..)
    -- , NormalErr(..)
    , ConfInt(..)
      -- ** Constructors
    -- , estimateNormErr
    , estimateFromInterval
    , estimateFromErr
      -- ** Accessors
    , confidenceInterval
    , Scale(..)
      -- * Other
    , Sample
    ) where

import Control.DeepSeq              (NFData(..))
import Data.Data                    (Data,Typeable)
import Data.Maybe                   (fromMaybe)
import GHC.Generics                 (Generic)

#if __GLASGOW_HASKELL__ == 704
import qualified Data.Vector.Generic
import qualified Data.Vector.Generic.Mutable
#endif

import Statistics.Internal
import Statistics.Types.Internal

----------------------------------------------------------------
-- Data type for confidence level
----------------------------------------------------------------

-- |
-- Confidence level. In context of confidence intervals it's
-- probability of said interval covering true value of measured
-- value. In context of statistical tests it's @1-α@ where α is
-- significance of test.
--
-- Since confidence level are usually close to 1 they are stored as
-- @1-CL@ internally. There are two smart constructors for @CL@:
-- 'mkCL' and 'mkCLFromSignificance' (and corresponding variant
-- returning @Maybe@). First creates @CL@ from confidence level and
-- second from @1 - CL@ or significance level.
--
-- >>> cl95
-- mkCLFromSignificance 0.05
--
-- Prior to 0.14 confidence levels were passed to function as plain
-- @Doubles@. Use 'mkCL' to convert them to @CL@.
newtype CL a = CL a
               deriving (Eq, Typeable, Data, Generic)

instance Show a => Show (CL a) where
  showsPrec n (CL p) = defaultShow1 "mkCLFromSignificance" p n
instance (Num a, Ord a, Read a) => Read (CL a) where
  readPrec = defaultReadPrecM1 "mkCLFromSignificance" mkCLFromSignificanceE

instance NFData   a => NFData   (CL a) where
  rnf (CL a) = rnf a

-- |
-- >>> cl95 > cl90
-- True
instance Ord a => Ord (CL a) where
  CL a <  CL b = a >  b
  CL a <= CL b = a >= b
  CL a >  CL b = a <  b
  CL a >= CL b = a <= b
  max (CL a) (CL b) = CL (min a b)
  min (CL a) (CL b) = CL (max a b)

-- | Create confidence level from probability β or probability
--   confidence interval contain true value of estimate. Will throw
--   exception if parameter is out of [0,1] range
--
-- >>> mkCL 0.95    -- same as cl95
-- mkCLFromSignificance 0.05
mkCL :: (Ord a, Num a) => a -> CL a
mkCL
  = fromMaybe (error "Statistics.Types.mkCL: probability is out if [0,1] range")
  . mkCLE

-- | Same as 'mkCL' but returns @Nothing@ instead of error if
--   parameter is out of [0,1] range
--
-- >>> mkCLE 0.95    -- same as cl95
-- Just (mkCLFromSignificance 0.05)
mkCLE :: (Ord a, Num a) => a -> Maybe (CL a)
mkCLE p
  | p >= 0 && p <= 1 = Just $ CL (1 - p)
  | otherwise        = Nothing

-- | Same as 'mkCLFromSignificance' but returns @Nothing@ instead of error if
--   parameter is out of [0,1] range
--
-- >>> mkCLFromSignificanceE 0.05    -- same as cl95
-- Just (mkCLFromSignificance 0.05)
mkCLFromSignificanceE :: (Ord a, Num a) => a -> Maybe (CL a)
mkCLFromSignificanceE p
  | p >= 0 && p <= 1 = Just $ CL p
  | otherwise        = Nothing

-- | Get confidence level. This function is subject to rounding
--   errors. If @1 - CL@ is needed use 'significanceLevel' instead
confidenceLevel :: (Num a) => CL a -> a
confidenceLevel (CL p) = 1 - p

-- | Get significance level.
significanceLevel :: CL a -> a
significanceLevel (CL p) = p

-- | 95% confidence level
cl95 :: Fractional a => CL a
cl95 = CL 0.05

----------------------------------------------------------------
-- Data type for p-value
----------------------------------------------------------------

-- | Newtype wrapper for p-value.
newtype PValue a = PValue a
               deriving (Eq,Ord, Typeable, Data, Generic)

instance Show a => Show (PValue a) where
  showsPrec n (PValue p) = defaultShow1 "mkPValue" p n
instance (Num a, Ord a, Read a) => Read (PValue a) where
  readPrec = defaultReadPrecM1 "mkPValue" mkPValueE

instance NFData a => NFData (PValue a) where
  rnf (PValue a) = rnf a

-- | Construct PValue. Returns @Nothing@ if argument is out of [0,1] range.
mkPValueE :: (Ord a, Num a) => a -> Maybe (PValue a)
mkPValueE p
  | p >= 0 && p <= 1 = Just $ PValue p
  | otherwise        = Nothing

----------------------------------------------------------------
-- Point estimates
----------------------------------------------------------------

-- |
-- A point estimate and its confidence interval. It's parametrized by
-- both error type @e@ and value type @a@. This module provides two
-- types of error: 'NormalErr' for normally distributed errors and
-- 'ConfInt' for error with normal distribution. See their
-- documentation for more details.
--
-- For example @144 ± 5@ (assuming normality) could be expressed as
--
-- > Estimate { estPoint = 144
-- >          , estError = NormalErr 5
-- >          }
--
-- Or if we want to express @144 + 6 - 4@ at CL95 we could write:
--
-- > Estimate { estPoint = 144
-- >          , estError = ConfInt
-- >                       { confIntLDX = 4
-- >                       , confIntUDX = 6
-- >                       , confIntCL  = cl95
-- >                       }
--
-- Prior to statistics 0.14 @Estimate@ data type used following definition:
--
-- > data Estimate = Estimate {
-- >      estPoint           :: {-# UNPACK #-} !Double
-- >    , estLowerBound      :: {-# UNPACK #-} !Double
-- >    , estUpperBound      :: {-# UNPACK #-} !Double
-- >    , estConfidenceLevel :: {-# UNPACK #-} !Double
-- >    }
--
-- Now type @Estimate ConfInt Double@ should be used instead. Function
-- 'estimateFromInterval' allow to easily construct estimate from same inputs.
data Estimate e a = Estimate
    { estPoint           :: !a
      -- ^ Point estimate.
    , estError           :: !(e a)
      -- ^ Confidence interval for estimate.
    } deriving (Eq, Read, Show, Generic
#if __GLASGOW_HASKELL__ >= 708
               , Typeable, Data
#endif
               )

instance (NFData   (e a), NFData   a) => NFData   (Estimate e a) where
    rnf (Estimate x dx) = rnf x `seq` rnf dx

-- | Confidence interval. It assumes that confidence interval forms
--   single interval and isn't set of disjoint intervals.
data ConfInt a = ConfInt
  { confIntLDX :: !a
    -- ^ Lower error estimate, or distance between point estimate and
    --   lower bound of confidence interval.
  , confIntUDX :: !a
    -- ^ Upper error estimate, or distance between point estimate and
    --   upper bound of confidence interval.
  , confIntCL  :: !(CL Double)
    -- ^ Confidence level corresponding to given confidence interval.
  }
  deriving (Read,Show,Eq,Typeable,Data,Generic)

instance NFData   a => NFData   (ConfInt a) where
    rnf (ConfInt x y _) = rnf x `seq` rnf y

----------------------------------------
-- Constructors

-- | Create estimate with asymmetric error.
estimateFromErr
  :: a                     -- ^ Central estimate
  -> (a,a)                 -- ^ Lower and upper errors. Both should be
                           --   positive but it's not checked.
  -> CL Double             -- ^ Confidence level for interval
  -> Estimate ConfInt a
estimateFromErr x (ldx,udx) cl = Estimate x (ConfInt ldx udx cl)

-- | Create estimate with asymmetric error.
estimateFromInterval
  :: Num a
  => a                     -- ^ Point estimate. Should lie within
                           --   interval but it's not checked.
  -> (a,a)                 -- ^ Lower and upper bounds of interval
  -> CL Double             -- ^ Confidence level for interval
  -> Estimate ConfInt a
estimateFromInterval x (lx,ux) cl
  = Estimate x (ConfInt (x-lx) (ux-x) cl)

----------------------------------------
-- Accessors

-- | Get confidence interval
confidenceInterval :: Num a => Estimate ConfInt a -> (a,a)
confidenceInterval (Estimate x (ConfInt ldx udx _))
  = (x - ldx, x + udx)

-- | Data types which could be multiplied by constant.
class Scale e where
  scale :: (Ord a, Num a) => a -> e a -> e a

instance Scale ConfInt where
  scale a (ConfInt l u cl) | a >= 0    = ConfInt  (a*l)  (a*u) cl
                           | otherwise = ConfInt (-a*u) (-a*l) cl

instance Scale e => Scale (Estimate e) where
  scale a (Estimate x dx) = Estimate (a*x) (scale a dx)

ext: hs

max_stars_repo_name: runeksvendsen/hs-gauge

max_issues_repo_name: runeksvendsen/hs-gauge

max_forks_repo_name: runeksvendsen/hs-gauge


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

-- Created September 7 2015

-- TODO | - UI and graphics settings
--        - Factor out logic (and settings, input mode, etc.) (should not be embedded in rendering functions)

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC Pragmas
--------------------------------------------------------------------------------------------------------------------------------------------
{-# LANGUAGE OverloadedStrings #-}

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module SpriteClip.Render (module SpriteClip.Render,
                          module Render,
                          module Palette,
                          module Southpaw.Picasso.RenderUtils,
                          module Cairo) where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Data.Complex
import Data.Maybe
import qualified Data.Set as S

import Text.Printf
import Control.Monad (forM_, when, liftM)
import Control.Lens

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

import           Southpaw.Picasso.Palette (alphapart, redpart, greenpart, bluepart)
import qualified Southpaw.Picasso.Palette     as Palette
import qualified Southpaw.Picasso.Render      as Render
import           Southpaw.Picasso.RenderUtils hiding (grid, closePath)

import qualified Southpaw.Cartesian.Plane.BoundingBox as BBox
import           Southpaw.Cartesian.Plane.Utilities

import           SpriteClip.Types
import           SpriteClip.Lenses
import qualified SpriteClip.Logic  as Logic

--------------------------------------------------------------------------------------------------------------------------------------------
-- Functions
--------------------------------------------------------------------------------------------------------------------------------------------
-- | Renders a frame
render :: AppState -> Cairo.Render ()
render appstate = do
  --
  canvasSize <- imageSurfaceSize sheet'

  -- Render.imageWithClip (const $ const $ Render.circle (appstate ^. mouse) 120) (50:+50) (appstate ^. sheet)
  Render.image (100:+100) (appstate ^. sheet)

  --
  -- Render.image (canvasSize*0.5) (appstate ^. sheet) -- Render sprite sheet
  -- Cairo.liftIO $ print (length $ take 5 sections)
  forM_ sections' (section mouse') -- Render sprite sections with markers
  alignments sections' mouse'      -- Render alignments
  axes mouse' canvasSize           -- Render cursor axes
  debugHUD appstate                -- Render debug overlay
  where
    sections' = (appstate ^. cutouts) ++ maybe [] (\click -> [BBox.fromCorners click mouse']) (appstate ^. pin)                -- TODO: Refactor
    mouse'    = Logic.applyOptionals [(gridsnap, dotmap $ Logic.nearestMultiple 20.0), (square, smallest)] (appstate ^. mouse) -- TODO: Rename (eg. cursor or pointer)
    sheet'    = appstate ^. sheet
    gridsnap  = S.member "shift_r"    $ appstate ^. keys --
    square    = S.member "control_r"  $ appstate ^. keys --
    smallest (re:+im) = min re im :+ min re im           --

-- |
-- TODO: Options
section ::  Complex Double -> BoundingBox Double -> Cairo.Render ()
section mouse bbox = do
  -- Outline
  Render.rectangle (_centre bbox) (_size bbox)
  choose Palette.magenta
  Cairo.stroke

  -- Markers
  sectionMarkers bbox mouse

-- | Renders the markers for a sheet section (respresented as a 'BoundingBox')
-- TODO: Options
-- TODO: Move out hover logic (?)
-- TODO: Refactor
sectionMarkers :: BoundingBox Double -> Complex Double -> Cairo.Render ()
sectionMarkers bbox mouse = forM_ markerBunches $ \ (radius, colour, offsets) ->
    forM_ offsets $ \pt -> do
      choose colour
      Render.circle (centre pt) (r radius pt)
      Cairo.fill
  where
    scale     = max 0.3 $ min 1.6 (let szx:+szy = _size bbox in min (abs szx) (abs szy)/120.0) -- Scale depending on section size
    focused r = Logic.within (scale*r) mouse . centre                                          -- Does the cursor lie on the marker?
    r r' pt   = scale*r'*(if focused r' pt then 1.4 else 1.0)                                  --
    (rectx:+recty)  = _size bbox
    centre (dx:+dy) = _centre bbox + ((dx*rectx/2):+(dy*recty/2))
    markerBunches   = [(12, Palette.orange & alphapart .~ 0.88, Logic.corners),
                       (10, Palette.peru   & alphapart .~ 0.88, Logic.midline),
                       ( 7, Palette.plum   & alphapart .~ 0.88, Logic.centre)]

-- | Renders the alignments for a sprite section (respresented as a 'BoundingBox')
-- TODO: Don't assume the point represents the mouse
alignments :: [BoundingBox Double] -> Complex Double -> Cairo.Render ()
alignments bboxes mouse@(mx:+my) = do

  -- TODO: Implement actual snapping
  -- TODO: Coalesce duplicate snaps

  forM_ (Logic.uniqueAlignments 0.0 6.0 bboxes mouse) (uncurry referenceline)

  -- TODO: Save and push state to prevent interference, rather than resetting manually
  Cairo.newPath        --
  Cairo.setDash [] 0.0 -- Disable dashes
  where
    crosshairs (fr:+om) = Cairo.setLineWidth 1.0 >> Cairo.setDash [] 0.0 >> choose Palette.black >> Render.crosshairs (fr:+om) (22:+22) >> Cairo.stroke
    cairosetup d        = choose (markTheSpot d) >> Cairo.setLineWidth 1 >> Cairo.setDash [8, 8] 0
    markTheSpot d       = let i = if d /= 0.0 then 0.4 else 0.0 in (i, i, i, 1.0) --

    referenceline (Vertical d)   (fr:+om) = cairosetup d >> Render.line (fr:+om) (fr:+my) >> Cairo.stroke >> crosshairs (fr:+om)
    referenceline (Horizontal d) (fr:+om) = cairosetup d >> Render.line (fr:+om) (mx:+om) >> Cairo.stroke >> crosshairs (fr:+om)

-- | Renders the X and Y axes meeting at the given origin point.
-- TODO: Options
-- TODO: Hide mouse, use axes instead
-- TODO: Use crosshairs with 'hollow middle' at the origin point (?)
axes :: Complex Double -> Complex Double -> Cairo.Render ()
axes (ox:+oy) (width:+height) = do
  Cairo.setLineWidth 1
  Cairo.setDash [5, 5] 14

  choose Palette.green
  Cairo.moveTo     0 oy
  Cairo.lineTo width oy
  Cairo.stroke

  choose Palette.red
  Cairo.moveTo ox      0
  Cairo.lineTo ox height
  Cairo.stroke

-- |
debugHUD :: AppState -> Cairo.Render ()
debugHUD appstate = do
  -- Debug info (HUD)
  Cairo.moveTo 20 20
  choose Palette.darkblue
  Cairo.setFontSize 16
  -- Cairo.fontOptionsSetAntilias Cairo.AntialiasDefault
  Cairo.selectFontFace ("Helvetica" :: String) Cairo.FontSlantNormal Cairo.FontWeightNormal
  Cairo.showText $ (printf "Mouse=(%.02f, %.02f)" mx my :: String)
  maybe pass (\(px:+py) -> Cairo.showText (printf " | Size=(%.02f, %.02f)" (abs $ px-mx) (abs $ py-my) :: String)) (appstate ^. pin)
  where
    (mx:+my) = appstate ^. mouse
    msize    = liftM (dotmap abs . subtract (mx:+my)) (appstate ^. pin)
    pass     = return ()

max_stars_repo_name: SwiftsNamesake/SpriteClip

max_issues_repo_name: SwiftsNamesake/SpriteClip

max_forks_repo_name: SwiftsNamesake/SpriteClip


text:
module Math.HiddenMarkovModel.Named (
   T(..),
   Discrete,
   Gaussian,
   fromModelAndNames,
   toCSV,
   fromCSV,
   ) where

import qualified Math.HiddenMarkovModel.Distribution as Distr
import qualified Math.HiddenMarkovModel.Private as HMM
import qualified Math.HiddenMarkovModel.CSV as HMMCSV
import Math.HiddenMarkovModel.Distribution (State(..))
import Math.HiddenMarkovModel.Utility (attachOnes)

import qualified Numeric.LinearAlgebra.Algorithms as Algo
import qualified Data.Packed.Vector as Vector

import qualified Text.CSV.Lazy.String as CSV
import Text.Printf (printf)

import qualified Control.Monad.Exception.Synchronous as ME
import qualified Control.Monad.Trans.State as MS
import Control.DeepSeq (NFData, rnf)
import Foreign.Storable (Storable)

import qualified Data.Map as Map
import qualified Data.List as List
import Data.Tuple.HT (swap)
import Data.Map (Map)

{- |
A Hidden Markov Model with names for each state.

Although 'nameFromStateMap' and 'stateFromNameMap' are exported
you must be careful to keep them consistent when you alter them.
-}
data T distr prob =
   Cons {
      model :: HMM.T distr prob,
      nameFromStateMap :: Map State String,
      stateFromNameMap :: Map String State
   }
   deriving (Show, Read)

type Discrete prob symbol = T (Distr.Discrete prob symbol) prob
type Gaussian a = T (Distr.Gaussian a) a

instance
   (NFData distr, NFData prob, Storable prob) =>
      NFData (T distr prob) where
   rnf hmm = rnf (model hmm, nameFromStateMap hmm, stateFromNameMap hmm)

fromModelAndNames :: HMM.T distr prob -> [String] -> T distr prob
fromModelAndNames md names =
   let m = Map.fromList $ zip [State 0 ..] names
   in  Cons {
          model = md,
          nameFromStateMap = m,
          stateFromNameMap = inverseMap m
       }

inverseMap :: Map State String -> Map String State
inverseMap =
   Map.fromListWith (error "duplicate label") .
   map swap . Map.toList

toCSV ::
   (Distr.CSV distr, Algo.Field prob, Show prob) =>
   T distr prob -> String
toCSV hmm =
   CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $
      Map.elems (nameFromStateMap hmm) : HMM.toCells (model hmm)

fromCSV ::
   (Distr.CSV distr, Algo.Field prob, Read prob) =>
   String -> ME.Exceptional String (T distr prob)
fromCSV =
   MS.evalStateT parseCSV . map HMMCSV.fixShortRow . CSV.parseCSV

parseCSV ::
   (Distr.CSV distr, Algo.Field prob, Read prob) =>
   HMMCSV.CSVParser (T distr prob)
parseCSV = do
   names <- HMMCSV.parseStringList =<< HMMCSV.getRow
   let duplicateNames =
         Map.keys $ Map.filter (> (1::Int)) $
         Map.fromListWith (+) $ attachOnes names
    in HMMCSV.assert (null duplicateNames) $
          "duplicate names: " ++ List.intercalate ", " duplicateNames
   md <- HMM.parseCSV
   let n = length names
       m = Vector.dim (HMM.initial md)
    in HMMCSV.assert (n == m) $
          printf "got %d state names for %d state" n m
   return $ fromModelAndNames md names

max_stars_repo_name: rybern/hmm-hmatrix

max_issues_repo_name: rybern/hmm-hmatrix

max_forks_repo_name: rybern/hmm-hmatrix


text:
{-# LANGUAGE LambdaCase #-}
module AI.NeuroCar where

import           AI.GeneticAlgorithm   as GA
import qualified AI.NeuralNetwork      as NN
import qualified Car                   as C
import           Control.Lens
import           Control.Monad.Loops
import           Control.Monad.State
import           Data.Foldable
import           Data.Functor.Identity
import           Data.Word             (Word32)
import           Game
import qualified Geometry              as G
import qualified Numeric.LinearAlgebra as LA
import           System.Random
import qualified Track                 as T
import           World

rayAngle :: Double
rayAngle = pi

rayCount :: Int
rayCount = 7

evolveCar :: Int -> Int -> Int -> Double -> Double -> Double -> Word32 ->
    C.CarParams -> T.Track Double -> [Population Double]
evolveCar seed generations popSize mutChance mutStrength time deltaTicks carParams track = do
    let gen          = mkStdGen seed
    let drawFunc     = return . const ()
    let inputFunc nn = Identity . getNetworkInput nn
    let timeFunc     = return deltaTicks
    let world        = initWorld carParams track time
    let evolveLoop nn = gameLoop drawFunc (inputFunc nn) timeFunc
    let runGame nn = runIdentity $ iterateUntilM (\w -> w^.gameState /= GameRunning) (evolveLoop nn) world
    let fitfunc nn = fromIntegral $ view score (runGame nn)
    let mutfunc = GA.mutate mutChance mutStrength
    let indGen = NN.newNetwork [3+rayCount, 15, 15, 2]
    let evoFunc = evolves generations popSize indGen fitfunc mutfunc
    evalState evoFunc gen

getNetworkInput :: NN.Network -> World -> [Input]
getNetworkInput nn w =
    let push a        = (state $ \xs -> ((),a:xs)) :: State [Input] ()
        rays          = C.shootRays rayAngle rayCount $ w^.car
        origin        = w^.car.C.position
        intersections = T.rayIntersection (w^.track) <$> rays
        distances     = (\case Nothing -> 0
                               Just v  -> G.segLength (origin, v)) <$> intersections
        velocity      = C.localVelocity $ w^.car
        rotation      = w^.car.C.rotation
        x             = rotation : toList velocity ++ distances
        y             = LA.toList $ NN.feedforward (LA.fromList x) nn
     in execState (do { when (head y < 1/3) (push GoBackward)
                      ; when (head y > 2/3) (push GoForward)
                      ; when (y !! 1 < 1/3) (push GoLeft)
                      ; when (y !! 1 > 2/3) (push GoRight) }) []

max_stars_repo_name: cornelius-sevald/neurocar

max_issues_repo_name: cornelius-sevald/neurocar

max_forks_repo_name: cornelius-sevald/neurocar


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies        #-}

import           Data.Bits
import           Data.Convertible
import           Data.List                            (maximumBy)
import           Data.Maybe                           (fromJust)
import           Data.Ord
import qualified Data.Vector.Storable                 as V
import           Data.Word
import           Foreign.ForeignPtr

import           Numeric.LinearAlgebra
import qualified Numeric.LinearAlgebra.Devel          as U
import qualified Numeric.LinearAlgebra.Static         as SA

import           Graphics.Gloss
import           Graphics.Gloss.Interface.IO.Game
import qualified Graphics.Gloss.Interface.IO.Interact as GI

import           Grenade
import           Grenade.Demos.MNIST                  hiding (runNet')
import           Grenade.Layers.Internal.Shrink

data MouseState = MouseDown | MouseUp

-- | all of the information needed to create the canvas to draw on
data World      = Canvas (S ('D2 224 224)) MouseState MNIST

-- | turns the matrix in the world structure into a picture. To do this, we use some pointer
--   magic to turn our matrix (which is represented internally as an array in memory) into
--   a bitmap, where each pixel is represented by four word8s for the R, G, B, A channels
--   (which we construct as a Word32)
renderCanvas :: World -> IO Picture
renderCanvas (Canvas (S2D !arr) _ _)
    = return $ bitmapOfForeignPtr 224 224 (BitmapFormat BottomToTop PxABGR) bitmapPtr False
  where
    -- | given a RealNum p, creates an RGBA word32 of the form (p, p, p, 255) - so it is essentially
    --   creating a greyscale pixel with no transparency.
    convColor :: RealNum -> Word32
    convColor p = let p'  =  convert p :: Word32
                      !w  =  unsafeShiftL p' 24
                         .|. unsafeShiftL p' 16
                         .|. unsafeShiftL p' 8
                         .|. 255
                  in w

    vec         = flatten $ SA.extract arr
    pixels      = V.map convColor vec
    (ptr, _, _) = U.unsafeToForeignPtr pixels
    -- | this last step is to turn the array of word32s into an array of word8s by
    --   interpreting each element as 4 word8s
    bitmapPtr   = castForeignPtr ptr

-- | updates the world with some event
handleInput :: GI.Event -> World -> IO World
-- if the mouse button is pressed down, update the mouse state
handleInput (GI.EventKey (GI.MouseButton GI.LeftButton) (GI.Down) _ _) (Canvas arr _ net)
  = return $ Canvas arr MouseDown net

-- if the mouse button is unpressed, update the mouse state and run a neural network inference
-- on the current canvas image by shrinking the image to the correct input size. Prints
-- the result to the console.
handleInput (GI.EventKey (GI.MouseButton GI.LeftButton) (GI.Up) _ _) (Canvas (S2D arr') _ net) =
  let extractedMatrix = SA.extract arr'
      shrunkenMatrix  = shrink_2d 224 224 28 28 extractedMatrix
      staticMatrix    = fromJust . fromStorableMatrix $ shrunkenMatrix
  in do
    putStrLn $ runNet' net staticMatrix
    return $ Canvas (S2D arr') MouseUp net

-- If the 'c' key is pressed and then let go of, wipe the canvas
handleInput (GI.EventKey (GI.Char 'c') (GI.Up) _ _) (Canvas (S2D !arr) mb net)
  = return $ Canvas cleanCanvas mb net
  where
    cleanCanvas = fromJust . fromStorableMatrix $ cleanMatrix
    cleanMatrix = U.mapMatrixWithIndex (const (const 255)) (SA.extract arr)

-- if moving the mouse while the left mouse button is held down, draw to the
-- canvas (note that we add 112, since the origin of a canvas is at the middle)
handleInput (GI.EventMotion (y, x)) (Canvas arr MouseDown net)
  = return $ Canvas (draw arr (xAsInt + 112) (yAsInt + 112)) MouseDown net
  where
    xAsInt = convert x
    yAsInt = convert y

handleInput _ c = return c

-- | draws a circle with squared radius 50 at the given coordinates
draw :: S ('D2 224 224) -> Int -> Int -> S ('D2 224 224)
draw (S2D arr) x y = fromJust . fromStorableMatrix $ m
  where
    m    = U.mapMatrixWithIndex f arr'
    arr' = SA.extract arr :: Matrix RealNum

    f :: (Int, Int) -> RealNum -> RealNum
    f (x', y') p = if (x - x') ^ (2 :: Int) + (y - y') ^ (2 :: Int) <= 50 then 0 else p

-- | runs a neural network inference and pretty prints the most probable label
--   with its probability
runNet' :: MNIST -> S ('D2 28 28) -> String
runNet' net m
  = let S1D ps = runNet net m
        (p, i) = (getProb . V.toList) (SA.extract ps)
    in  "This number is " ++ show i ++ " with probability " ++ show (p * 100) ++ "%"
  where
    getProb :: [RealNum] -> (RealNum, Int)
    getProb xs = maximumBy (comparing fst) (zip xs [0..])

main :: IO ()
main = do
    mnistPath <- getPathForNetwork MNIST
    net <- (loadSerializedNetwork mnistPath :: IO MNIST)
    putStrLn "Successfully loaded model"

    let initialCanvas = Canvas initialMat MouseUp net

    playIO window backgroundColor 30 initialCanvas renderCanvas handleInput emptyStep
  where
    window            = InWindow "Draw here!" (224, 224) (100, 100)
    backgroundColor   = makeColor 255 255 255 0
    initialMat        = fromJust . fromStorable . V.fromList $ replicate (224 * 224) 255
    emptyStep _ world = return world

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-|
Module: MachineLearning.NeuralNetwork.SoftmaxLoss
Description: Softmax Loss.
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX

Softmax Loss.
-}

module MachineLearning.NeuralNetwork.SoftmaxLoss
(
  scores
  , gradient
  , loss
)

where

import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA
import MachineLearning.Types (R, Matrix)
import MachineLearning.Utils (sumByRows, reduceByRows)

scores x = x - reduceByRows V.maximum x

gradient scores y =
  let sum_probs = sumByRows $ exp scores
      probs = (exp scores) / sum_probs
  in probs - y

-- Softmax Loss function
loss :: Matrix -> Matrix -> R
loss scores y = loss / m
  where m = fromIntegral $ LA.rows scores
        sum_probs = sumByRows $ exp scores
        t = sumByRows $ scores * y
        loss = LA.sumElements $ (log sum_probs) - t

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Polynomial (multiply, multiplyWithIntCoeffs, Poly(Poly), Coefficients) where

import Data.Complex
import DFT (dft, idft)
import Data.List (take, zip, repeat)

type Coefficients = [Complex Float]
data Poly = Poly Coefficients deriving (Eq)

instance Show Poly where
  show (Poly p) = show' p 0
    where show' [] degree = ""
          show' (x:xs) 0 = "(" ++ show (round (realPart x)) ++ " + " ++ ("(" ++ show (round (imagPart x)) ++ ")") ++ "i" ++ ")" ++ " + " ++ show' xs 1
          show' [x] degree = "(" ++ show (round (realPart x)) ++ " + " ++ ("(" ++ show (round (imagPart x)) ++ ")") ++ "i" ++ ")" ++ "x^" ++ show degree
          show' (x:xs) degree = "(" ++ show (round (realPart x)) ++ " + " ++ ("(" ++ show (round (imagPart x)) ++ ")") ++ "i" ++ ")" ++ "x^" ++ show degree ++ " + " ++ show' xs (degree + 1)

multiplyWithIntCoeffs :: Poly -> Poly -> [Int]
multiplyWithIntCoeffs p q = map round $ map (realPart) $ coeffs
  where (Poly coeffs) = multiply p q

multiply :: Poly -> Poly -> Poly
multiply (Poly p) (Poly q) = Poly (trimPoly $ idft (dft p' |*| dft q'))
  where (p', q') = pad p q
        trimPoly = take (length p + length q - 1)

pad :: Coefficients -> Coefficients -> (Coefficients, Coefficients)
pad p q = (p', q')
  where maxDegree = closestPowerOfTwo (length p + length q + 2)
        p' = p ++ (take (maxDegree - (length p)) (repeat 0))
        q' = q ++ (take (maxDegree - (length q)) (repeat 0))

(|*|) :: Coefficients -> Coefficients -> Coefficients
(|*|) = hadamard

hadamard :: Coefficients -> Coefficients -> Coefficients
hadamard p q = [ a * b | (a, b) <- zip p q ]

closestPowerOfTwo :: Int -> Int
closestPowerOfTwo x = 2 ^ ceiling (log (fromIntegral x) / log 2)

max_stars_repo_name: gboduljak/dft

max_issues_repo_name: gboduljak/dft

max_forks_repo_name: gboduljak/dft


text:
module Statistics.Iteratee.Tests

where

import Data.Iteratee as I
import Statistics.Iteratee as Si
import Statistics.Sample as St
import Test.Framework
import Test.Framework.Providers.QuickCheck2 (testProperty)

import Control.Monad.Identity
import qualified Data.Vector.Unboxed as V

tests =
  [ testGroup "Sample" $ map mkUnProp uns
  ]

-- Don't want to drag in more dependencies for just this.
class ApproxEq a where
    approxEq :: Double -> a -> a -> Bool

instance ApproxEq Double where
    approxEq tol d1 d2 = m == 0.0 || d/m < tol where
      m = max (abs d1) (abs d2)
      d = abs (d1 - d2)

instance (ApproxEq a, ApproxEq b) => ApproxEq (a,b) where
    approxEq tol (l1,r1) (l2,r2) = approxEq tol l1 l2 && approxEq tol r1 r2

infix 4 ===
(===) :: (ApproxEq a) => a -> a -> Bool
(===) = approxEq {-pretty equal-}1.0e-10

unsProp :: (ApproxEq a)
        => (V.Vector Double -> a)
        -> (Iteratee [Double] Identity a)
        -> [Double]
        -> Bool
unsProp vec iter xs = if null xs then True
    else vec (V.fromList xs) === (runIdentity $ run =<< enumPure1Chunk xs iter)
    -- we're using Eq for doubles, which is always a bad idea...

    -- also not checking empty vectors, because in some cases (range,
    -- harmonicMean) we get NaN's or other funky values.

mkUnProp (lbl, st, si) = testProperty lbl $ unsProp st si

-- unary properties
uns =
  [ ("mean", St.mean, Si.mean)
  , ("range", St.range, Si.range)
  , ("harmonic_mean", St.harmonicMean, Si.harmonicMean)
  , ("variance", St.fastVariance, Si.variance)
  , ("std_dev", St.fastStdDev,   Si.stdDev)
  ]

max_stars_repo_name: JohnLato/iter-stats

max_issues_repo_name: JohnLato/iter-stats

max_forks_repo_name: JohnLato/iter-stats


text:
module Statistics.Information.Mixed.Entropy where

import Data.Matrix
import qualified Statistics.Information.Continuous.Entropy as CE
import qualified Statistics.Information.Discrete.Entropy as DE
import Statistics.Information.Mixed.MutualInfo
import Statistics.Information.Utils.Random
import System.Random

centropydc :: (RandomGen g, Eq a) => g -> Int -> Int -> Matrix a ->
              Matrix Double -> Double
centropydc g k base xs ys = DE.entropy base xs - midc g k base xs ys

centropycd :: (RandomGen g, Eq a) => g -> Int -> Int -> Matrix Double ->
              Matrix a -> Double
centropycd g k base xs ys = CE.entropy g1 k base xs - micd g2 k base xs ys where
  (g1, g2) = split g

entropy :: (RandomGen g, Eq a) => g -> Int -> Int -> Matrix a ->
           Matrix Double -> Double
entropy g k base xs ys = hx + hy - ixy
  where
    (g1, g2) = split g
    hx = DE.entropy base xs
    hy = CE.entropy g1 k base ys
    ixy = midc g2 k base xs ys

max_stars_repo_name: eligottlieb/Shannon

max_issues_repo_name: eligottlieb/Shannon

max_forks_repo_name: eligottlieb/Shannon


text:
import Data.List(transpose)
import Numeric.LinearAlgebra.Algorithms

data Matrix a = Matrix { width :: Int, height :: Int, array :: [a] } deriving (Show, Eq)

mkMatrix :: (Num a) => Int -> Int -> [a] -> Matrix a
mkMatrix w h a = Matrix w h (take (w*h) (a ++ repeat 0))

--add :: (Num a, Num b, Num c) => Matrix a -> Matrix b -> Maybe (Matrix c)
add (Matrix w h a) (Matrix w' h' a')
  | w == w' && h == h' = Just $ Matrix w h (zipWith (+) a a')
  | otherwise = Nothing

multiply (Matrix w h a) (Matrix w' h' a')
  | w == h' = undefined --normal
  | otherwise = Nothing

groupN n xs@(_:_) = (take n xs) : groupN n (drop n xs)
groupN n [] = []

transpose' :: [[a]] -> [[a]]
transpose' xs@(x:_)
  | null x = []
  | otherwise = heads : transpose' tails
    where heads = fmap head xs
          tails = fmap tail xs


text:
{-# LANGUAGE
    MultiWayIf,
    FlexibleInstances,
    NumDecimals #-}

module Jac where

import Data.Matrix (Matrix(..), fromLists, fromList, toList)
import Data.Complex

import Point
import Data (F)

type FPD = Point -> Double

-----------------------------------------------------------------
-- производная.
diff :: (Floating a, Eq a) => F a -> a -> a
diff f x = (f (x - d/2) - f (x + d/2))/d
  where
    d = if
        | x == 0    ->   1e-9
        | otherwise -> x*1e-5

-- якобиан.
jac :: F Point -> Point -> Matrix Double
jac f p@(Point x' y') = fromLists [
    [diff (\a -> f1 p{x = a}) x', diff (\a -> f1 p{y = a}) y'],
    [diff (\a -> f2 p{x = a}) x', diff (\a -> f2 p{y = a}) y']]
  where
    f1 = x.f ; f2 = y.f

jac' :: F Point -> Point -> Double
jac' f p@(Point x' y') = log $ abs $
    diff (\a -> f1 p{x = a}) x'*diff (\a -> f2 p{y = a}) y'
  - diff (\a -> f1 p{y = a}) y'*diff (\a -> f2 p{x = a}) x'
  where
    f1 = x.f ; f2 = y.f

-- единичный вектор.
e :: Double -> Matrix Double
e f = fromList 2 1 [realPart x, imagPart x]
  where x = cis f

-- см. стр. 102. в книге 1.
formA :: F Point -> Point -> Double -> Double
formA f p k = log.abs.toList $ jac f p * e k
    where abs [x, y] = magnitude $ x :+ y

-----------------------------------------------------------------
-- преобразование вектора в точку.
instance Points (Matrix Double) where
    toPoint = toPoint.toList
    fromPoint = fromList 2 1.fromPoint

-----------------------------------------------------------------
--                   Справочная литература                     --
-----------------------------------------------------------------
-- 1. Г.С. Осипенко, Н.Б.Ампилова:
--      ВВЕДЕНИЕ В СИМВОЛИЧЕСКИЙ АНАЛИЗ ДИНАМИЧЕСКИХ СИСТЕМ

max_stars_repo_name: vojiranto/dsa

max_issues_repo_name: vojiranto/dsa

max_forks_repo_name: vojiranto/dsa


text:
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.Lognormal
-- Copyright : (c) 2020 Ximin Luo
-- License   : BSD3
--
-- Maintainer  : infinity0@pwned.gg
-- Stability   : experimental
-- Portability : portable
--
-- The Weibull distribution.  This is a continuous probability
-- distribution that describes the occurrence of a single event whose
-- probability changes over time, controlled by the shape parameter.

module Statistics.Distribution.Weibull
    (
      WeibullDistribution
      -- * Constructors
    , weibullDistr
    , weibullDistrErr
    , weibullStandard
    , weibullDistrApproxMeanStddevErr
    ) where

import Control.Applicative
import Data.Aeson            (FromJSON(..), ToJSON, Value(..), (.:))
import Data.Binary           (Binary(..))
import Data.Data             (Data, Typeable)
import GHC.Generics          (Generic)
import Numeric.MathFunctions.Constants (m_eulerMascheroni)
import Numeric.SpecFunctions (expm1, log1p, logGamma)
import qualified Data.Vector.Generic as G

import qualified Statistics.Distribution as D
import qualified Statistics.Sample as S
import Statistics.Internal

-- | The Weibull distribution.
data WeibullDistribution = WD {
      wdShape  :: {-# UNPACK #-} !Double
    , wdLambda :: {-# UNPACK #-} !Double
    } deriving (Eq, Typeable, Data, Generic)

instance Show WeibullDistribution where
  showsPrec i (WD k l) = defaultShow2 "weibullDistr" k l i
instance Read WeibullDistribution where
  readPrec = defaultReadPrecM2 "weibullDistr" $
    (either (const Nothing) Just .) . weibullDistrErr

instance ToJSON WeibullDistribution
instance FromJSON WeibullDistribution where
  parseJSON (Object v) = do
    k <- v .: "wdShape"
    l <- v .: "wdLambda"
    either fail return $ weibullDistrErr k l
  parseJSON _ = empty

instance Binary WeibullDistribution where
  put (WD k l) = put k >> put l
  get = do
    k <- get
    l <- get
    either fail return $ weibullDistrErr k l

instance D.Distribution WeibullDistribution where
  cumulative      = cumulative
  complCumulative = complCumulative

instance D.ContDistr WeibullDistribution where
  logDensity    = logDensity
  quantile      = quantile
  complQuantile = complQuantile

instance D.MaybeMean WeibullDistribution where
  maybeMean = Just . D.mean

instance D.Mean WeibullDistribution where
  mean (WD k l) = l * exp (logGamma (1 + 1 / k))

instance D.MaybeVariance WeibullDistribution where
  maybeStdDev   = Just . D.stdDev
  maybeVariance = Just . D.variance

instance D.Variance WeibullDistribution where
  variance (WD k l) = l * l * (exp (logGamma (1 + 2 * invk)) - q * q)
   where
    invk = 1 / k
    q    = exp (logGamma (1 + invk))

instance D.Entropy WeibullDistribution where
  entropy (WD k l) = m_eulerMascheroni * (1 - 1 / k) + log (l / k) + 1

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

instance D.ContGen WeibullDistribution where
  genContVar d = D.genContinuous d

-- | Standard Weibull distribution with scale factor (lambda) 1.
weibullStandard :: Double -> WeibullDistribution
weibullStandard k = weibullDistr k 1.0

-- | Create Weibull distribution from parameters.
--
-- If the shape (first) parameter is @1.0@, the distribution is equivalent to a
-- 'Statistics.Distribution.Exponential.ExponentialDistribution' with parameter
-- @1 / lambda@ the scale (second) parameter.
weibullDistr
  :: Double            -- ^ Shape
  -> Double            -- ^ Lambda (scale)
  -> WeibullDistribution
weibullDistr k l = either error id $ weibullDistrErr k l

-- | Create Weibull distribution from parameters.
--
-- If the shape (first) parameter is @1.0@, the distribution is equivalent to a
-- 'Statistics.Distribution.Exponential.ExponentialDistribution' with parameter
-- @1 / lambda@ the scale (second) parameter.
weibullDistrErr
  :: Double            -- ^ Shape
  -> Double            -- ^ Lambda (scale)
  -> Either String WeibullDistribution
weibullDistrErr k l | k <= 0     = Left $ errMsg k l
                    | l <= 0     = Left $ errMsg k l
                    | otherwise = Right $ WD k l

errMsg :: Double -> Double -> String
errMsg k l =
  "Statistics.Distribution.Weibull.weibullDistr: both shape and lambda must be positive. Got shape "
    ++ show k
    ++ " and lambda "
    ++ show l

-- | Create Weibull distribution from mean and standard deviation.
--
-- The algorithm is from "Methods for Estimating Wind Speed Frequency
-- Distributions", C. G. Justus, W. R. Hargreaves, A. Mikhail, D. Graber, 1977.
-- Given the identity:
--
-- \[
-- (\frac{\sigma}{\mu})^2 = \frac{\Gamma(1+2/k)}{\Gamma(1+1/k)^2} - 1
-- \]
--
-- \(k\) can be approximated by
--
-- \[
-- k \approx (\frac{\sigma}{\mu})^{-1.086}
-- \]
--
-- \(\lambda\) is then calculated straightforwardly via the identity
--
-- \[
-- \lambda = \frac{\mu}{\Gamma(1+1/k)}
-- \]
--
-- Numerically speaking, the approximation for \(k\) is accurate only within a
-- certain range. We arbitrarily pick the range \(0.033 \le \frac{\sigma}{\mu} \le 1.45\)
-- where it is good to ~6%, and will refuse to create a distribution outside of
-- this range. The paper does not cover these details but it is straightforward
-- to check them numerically.
weibullDistrApproxMeanStddevErr
  :: Double            -- ^ Mean
  -> Double            -- ^ Stddev
  -> Either String WeibullDistribution
weibullDistrApproxMeanStddevErr m s = if r > 1.45 || r < 0.033
    then Left msg
    else weibullDistrErr k l
  where r = s / m
        k = (s / m) ** (-1.086)
        l = m / exp (logGamma (1 + 1/k))
        msg = "Statistics.Distribution.Weibull.weibullDistr: stddev-mean ratio "
          ++ "outside approximation accuracy range [0.033, 1.45]. Got "
          ++ "stddev " ++ show s ++ " and mean " ++ show m

-- | Uses an approximation based on the mean and standard deviation in
--   'weibullDistrEstMeanStddevErr', with standard deviation estimated
--   using maximum likelihood method (unbiased estimation).
--
--   Returns @Nothing@ if sample contains less than one element or
--   variance is zero (all elements are equal), or if the estimated mean
--   and standard-deviation lies outside the range for which the
--   approximation is accurate.
instance D.FromSample WeibullDistribution Double where
  fromSample xs
    | G.length xs <= 1 = Nothing
    | v == 0           = Nothing
    | otherwise        = either (const Nothing) Just $
      weibullDistrApproxMeanStddevErr m (sqrt v)
    where
      (m,v) = S.meanVarianceUnb xs

logDensity :: WeibullDistribution -> Double -> Double
logDensity (WD k l) x
  | x < 0     = 0
  | otherwise = log k + (k - 1) * log x - k * log l - (x / l) ** k

cumulative :: WeibullDistribution -> Double -> Double
cumulative (WD k l) x | x < 0     = 0
                      | otherwise = -expm1 (-(x / l) ** k)

complCumulative :: WeibullDistribution -> Double -> Double
complCumulative (WD k l) x | x < 0     = 1
                           | otherwise = exp (-(x / l) ** k)

quantile :: WeibullDistribution -> Double -> Double
quantile (WD k l) p
  | p == 0         = 0
  | p == 1         = inf
  | p > 0 && p < 1 = l * (-log1p (-p)) ** (1 / k)
  | otherwise      =
    error $ "Statistics.Distribution.Weibull.quantile: p must be in [0,1] range. Got: " ++ show p
  where inf = 1 / 0

complQuantile :: WeibullDistribution -> Double -> Double
complQuantile (WD k l) q
  | q == 0         = inf
  | q == 1         = 0
  | q > 0 && q < 1 = l * (-log q) ** (1 / k)
  | otherwise      =
    error $ "Statistics.Distribution.Weibull.complQuantile: q must be in [0,1] range. Got: " ++ show q
  where inf = 1 / 0

max_stars_repo_name: vaerksted/statistics

max_issues_repo_name: vaerksted/statistics

max_forks_repo_name: vaerksted/statistics


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Math.Tensor.Internal.LinearAlgebra
-- Copyright   :  (c) 2019 Tobias Reinhart and Nils Alex
-- License     :  MIT
-- Maintainer  :  tobi.reinhart@fau.de, nils.alex@fau.de
--
-- Gaussian elimination algorithm based on hmatrix.
-----------------------------------------------------------------------------

{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleContexts #-}

module Math.Tensor.Internal.LinearAlgebra (
-- * Gaussian Elimination
gaussianST, gaussianFFST,
gaussian, gaussianFF,
rrefST,
rref,
-- * Linearly Independent Columns
independentColumns, independentColumnsFF, independentColumnsRREF,
independentColumnsVerifiedFF,
independentColumnsMat, independentColumnsMatFF, independentColumnsMatRREF,
-- * Pivots
pivotsU, pivotsUFF,
findPivotMax, findPivotMaxFF, findRowPivot,
isref, isrref, isrref', verify)

where

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data
import Numeric.LinearAlgebra.Devel

import Debug.Trace (trace)

import Data.List (maximumBy)

import Control.Monad
import Control.Monad.ST

-- | Returns the pivot columns of an upper triangular matrix.
--
-- @
-- λ let mat = (3 >< 4) [1, 0, 2, -3, 0, 0, 1, 0, 0, 0, 0, 0]
-- λ mat
-- (3><4)
--  [ 1.0, 0.0, 2.0, -3.0
--  , 0.0, 0.0, 1.0,  0.0
--  , 0.0, 0.0, 0.0,  0.0 ]
-- λ pivotsU mat
-- [0,2]
-- @
--

pivotsU :: Matrix Double -> [Int]
pivotsU mat = go (0,0)
  where
    go (i,j)
      = case findPivot mat e (i,j) of
          Nothing       -> []
          Just (i', j') -> j' : go (i'+1, j'+1)
    maxAbs = maximum $ map (maximum . map abs) $ toLists mat
    e = eps * maxAbs

pivotsUFF :: Matrix Z -> [Int]
pivotsUFF mat = go (0,0)
  where
    go (i,j)
      = case findPivotFF mat (i,j) of
          Nothing       -> []
          Just (i', j') -> j' : go (i'+1, j'+1)

eps :: Double
eps = 1e-12

-- find next pivot in upper triangular matrix

findPivotFF :: Matrix Z -> (Int, Int) -> Maybe (Int, Int)
findPivotFF mat (i, j)
    | n == j = Nothing
    | m == i = Nothing
    | otherwise = case nonZeros of
                    []          -> if n == j+1
                                   then Nothing
                                   else findPivotFF mat (i, j+1)
                    (pi, pj):_  -> Just (pi, pj+j)
    where
        m = rows mat
        n = cols mat
        col = mat ¿ [j]
        nonZeros = filter (\(i', _) -> i' >= i) $ find (/= 0) col

findPivot :: Matrix Double -> Double -> (Int, Int) -> Maybe (Int, Int)
findPivot mat e (i, j)
    | n == j = Nothing
    | m == i = Nothing
    | otherwise = case nonZeros of
                    []          -> if n == j+1
                                   then Nothing
                                   else findPivot mat e (i, j+1)
                    (pi, pj):_  -> Just (pi, pj+j)
    where
        m = rows mat
        n = cols mat
        col = mat ¿ [j]
        nonZeros = filter (\(i', _) -> i' >= i) $ find (not . (< e) . abs) col

-- | Find pivot element below position (i, j) with greatest absolute value in the ST monad.

findPivotMaxFF :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe (Int, Int))
findPivotMaxFF m n i j mat
    | n == j = return Nothing
    | m == i = return Nothing
    | otherwise =
        do
          col      <- mapM (\i' -> do
                                    x <- readMatrix mat i' j
                                    return (i', x))
                      [i..m-1]
          let nonZeros = filter ((/= 0) . snd) col
          case nonZeros of
            []        -> if n == j+1
                         then return Nothing
                         else findPivotMaxFF m n i (j+1) mat
            (pi,_):_  -> return $ Just (pi, j)

findPivotMax :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s (Maybe (Int, Int))
findPivotMax m n i j mat
    | n == j = return Nothing
    | m == i = return Nothing
    | otherwise =
        do
          col      <- mapM (\i' -> do
                                    x <- readMatrix mat i' j
                                    return (i', abs x))
                      [i..m-1]
          let nonZeros = filter (not . (<eps) . abs . snd) col
          let (pi, _) = maximumBy (\(_, x) (_, y) -> x `compare` y) nonZeros
          case nonZeros of
            [] -> if n == j+1
                  then return Nothing
                  else findPivotMax m n i (j+1) mat
            _  -> return $ Just (pi, j)

-- gaussian elimination of sub matrix below position (i, j)

findRowPivot :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe Int)
findRowPivot m n i j mat
    | j + 1 > n       = error "out of bounds" -- return Nothing
    | i + 1 > min m n = error "out of bounds" -- return Nothing
    | otherwise =
        do
         row <- mapM (\j' -> do
                              x <- readMatrix mat i j'
                              return (j', x))
                [0 .. j]
         let nonZeros = filter ((/=0) . snd) row
         case nonZeros of
           []        -> return Nothing
           (pj, _):_ -> return $ Just pj

backwardFF' :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
backwardFF' m n i j mat
      | i == 0 = return ()
      | otherwise = do
    mat' <- freezeMatrix mat
    iPivot' <- findRowPivot m n i j mat
    case iPivot' of
        Nothing -> backwardFF' m n (i-1) j mat
        Just pj -> do
                    pv <- readMatrix mat i pj
                    mapM_ (reduce pv pj) [0 .. i-1]
                    backwardFF' m n (i-1) (pj-1) mat
  where
    reduce pv pj r = do
                      Just pr <- findRowPivot m n r pj mat
                      -- prv <- readMatrix mat r pr
                      pjv <- readMatrix mat r pj
                      if pjv == 0
                        then return ()
                        else
                         let op1 = SCAL pv (Row r) (FromCol pr)
                             op2 = AXPY (-pjv) i r (FromCol pj)
                         in do
                             rowOper op1 mat
                             rowOper op2 mat
                             g <- foldM (\acc c -> fmap (fromIntegral . gcd acc . fromIntegral) $ readMatrix mat r c) 0 [pr .. n-1]
                             if g == 0
                               then return()
                               else mapM_ (\c -> modifyMatrix mat r c (`quot` g)) [pr .. n-1]

gaussianFF' :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
gaussianFF' m n i j mat = do
    iPivot' <- findPivotMaxFF m n i j mat
    case iPivot' of
        Nothing     -> return ()
        Just (r, p) -> do
                          rowOper (SWAP i r (FromCol j)) mat
                          pv <- readMatrix mat i p
                          mapM_ (reduce pv p) [i+1 .. m-1]
                          gaussianFF' m n (i+1) (p+1) mat
  where
    reduce pv p r = do
                      rv <- readMatrix mat r p
                      if rv == 0
                        then return ()
                        else
                         let op1 = SCAL pv (Row r) (FromCol p)
                             op2 = AXPY (-rv) i r (FromCol p)
                         in do
                             rowOper op1 mat
                             rowOper op2 mat
                             g <- foldM (\acc c -> fmap (fromIntegral . gcd acc . fromIntegral) $ readMatrix mat r c) 0 [p .. n-1]
                             if g == 0
                               then return()
                               else mapM_ (\c -> modifyMatrix mat r c (`quot` g)) [p .. n-1]

gaussian' :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s ()
gaussian' m n i j mat = do
    iPivot' <- findPivotMax m n i j mat
    case iPivot' of
        Nothing     -> return ()
        Just (r, p) -> do
                          rowOper (SWAP i r (FromCol j)) mat
                          pv <- readMatrix mat i p
                          mapM_ (reduce pv p) [i+1 .. m-1]
                          gaussian' m n (i+1) (p+1) mat
  where
    reduce pv p r = do
                      rv <- readMatrix mat r p
                      if abs rv < eps
                        then return ()
                        else
                         let frac = -rv / pv
                             op = AXPY frac i r (FromCol p)
                         in do
                             rowOper op mat
                             mapM_ (\j' -> modifyMatrix mat r j' (\x -> if abs x < eps then 0 else x)) [p..n-1]

-- | Gaussian elimination perfomed in-place in the @'ST'@ monad.

gaussianFFST :: Int -> Int -> STMatrix s Z -> ST s ()
gaussianFFST m n = gaussianFF' m n 0 0

gaussianST :: Int -> Int -> STMatrix s Double -> ST s ()
gaussianST m n = gaussian' m n 0 0

rrefST :: Int -> Int -> STMatrix s Z -> ST s ()
rrefST m n mat = do
                    gaussianFF' m n 0 0 mat
                    backwardFF' m n (r'-1) (n-1) mat
    where
        r' = min m n

-- | Gaussian elimination as pure function. Involves a copy of the input matrix.
--
-- @
-- λ let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1]
-- λ mat
-- (3><4)
--  [ 1.0, 1.0, -2.0,  0.0
--  , 0.0, 2.0, -6.0, -4.0
--  , 3.0, 0.0,  3.0,  1.0 ]
-- λ gaussian mat
-- (3><4)
--  [ 3.0, 0.0,  3.0,                1.0
--  , 0.0, 2.0, -6.0,               -4.0
--  , 0.0, 0.0,  0.0, 1.6666666666666667 ]
-- @
--

isref :: (Element a, Num a, Ord a, Container Vector a) => Matrix a -> Bool
isref mat = case pivot of
              []      -> True
              (r,p):_ -> if r > 0
                         then False
                         else let leftMat  = mat ?? (Drop 1, Take (p+1))
                                  rightMat = mat ?? (Drop 1, Drop (p+1))
                                  leftZero = null $ find (/=0) leftMat
                                  rightRef = isref rightMat
                              in leftZero && rightRef
    where
        pivot = find (/=0) mat

isrref' :: (Element a, Num a, Ord a, Container Vector a) => Int -> Matrix a -> Bool
isrref' r mat = case pivot of
              []       -> True
              (r',p):_ -> if r' > 0
                          then False
                          else let leftMat  = subMat ?? (Drop 1, Take (p+1))
                                   col      = mat ¿ [p]
                                   colSingleton = case find (/=0) col of
                                                    [_] -> True
                                                    _   -> False
                                   leftZero = null $ find (/=0) leftMat
                                   nextRref = isrref' (r+1) mat
                               in leftZero && colSingleton && nextRref
    where
        subMat = mat ?? (Drop r, All)
        pivot  = find (/=0) subMat

isrref :: (Element a, Num a, Ord a, Container Vector a) => Matrix a -> Bool
isrref = isrref' 0

rref :: Matrix Z -> Matrix Z
rref mat = runST $ do
    matST <- thawMatrix mat
    rrefST m n matST
    freezeMatrix matST
  where
    m = rows mat
    n = cols mat

gaussianFF :: Matrix Z -> Matrix Z
gaussianFF mat = runST $ do
    matST <- thawMatrix mat
    gaussianFFST m n matST
    freezeMatrix matST
  where
    m = rows mat
    n = cols mat

gaussian :: Matrix Double -> Matrix Double
gaussian mat = runST $ do
    matST <- thawMatrix mat
    gaussianST m n matST
    freezeMatrix matST
  where
    m = rows mat
    n = cols mat

-- | Returns the indices of a maximal linearly independent subset of the columns
--   in the matrix.
--
-- @
-- λ let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1]
-- λ mat
-- (3><4)
--  [ 1.0, 1.0, -2.0,  0.0
--  , 0.0, 2.0, -6.0, -4.0
--  , 3.0, 0.0,  3.0,  1.0 ]
-- λ independentColumns mat
-- [0,1,3]
-- @
--

independentColumnsRREF :: Matrix Z -> [Int]
independentColumnsRREF mat = pivotsUFF mat'
    where
        mat' = rref mat

independentColumnsFF :: Matrix Z -> [Int]
independentColumnsFF mat = pivotsUFF mat'
    where
        mat' = gaussianFF mat

independentColumnsVerifiedFF :: Matrix Z -> [Int]
independentColumnsVerifiedFF mat
        | isref mat' && verify mat mat'
                     = pivotsUFF mat'
        | otherwise  = error "could not verify row echelon form"
    where
        mat' = gaussianFF mat

independentColumns :: Matrix Double -> [Int]
independentColumns mat = pivotsU mat'
    where
        mat' = gaussian mat

-- | Returns a sub matrix containing a maximal linearly independent subset of
--   the columns in the matrix.
--
-- @
-- λ let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1]
-- λ mat
-- (3><4)
--  [ 1.0, 1.0, -2.0,  0.0
--  , 0.0, 2.0, -6.0, -4.0
--  , 3.0, 0.0,  3.0,  1.0 ]
-- λ independentColumnsMat mat
-- (3><3)
--  [ 1.0, 1.0,  0.0
--  , 0.0, 2.0, -4.0
--  , 3.0, 0.0,  1.0 ]
-- @
--

verify :: Matrix Z -> Matrix Z -> Bool
verify mat ref = rank1 == rank2 && rank1 == rank3
    where
        matD = fromZ mat :: Matrix Double
        refD = fromZ ref :: Matrix Double
        rank1 = rank matD
        rank2 = rank refD
        rank3 = rank $ matD === refD

independentColumnsMatRREF :: Matrix Z -> Matrix Z
independentColumnsMatRREF mat =
  case independentColumnsRREF mat of
    [] -> (rows mat >< 1) $ repeat 0
    cs -> mat ¿ cs

independentColumnsMatFF :: Matrix Z -> Matrix Z
independentColumnsMatFF mat =
  case independentColumnsFF mat of
    [] -> (rows mat >< 1) $ repeat 0
    cs -> mat ¿ cs

independentColumnsMat :: Matrix Double -> Matrix Double
independentColumnsMat mat =
  case independentColumns mat of
    [] -> (rows mat >< 1) $ repeat 0
    cs -> mat ¿ cs


text:
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Commons where

import Control.Applicative (liftA2)
import Control.Monad (foldM, forM)
import Data.Array
import Data.Complex
import Data.Function.HT (nest)
import qualified Data.IntMap.Strict as IM
import Data.List (intercalate, sort)
import Data.Map.Strict (Map)
import qualified Data.Map.Strict as Map
import Data.Maybe (catMaybes, fromJust, mapMaybe)
import Data.Set (Set, fromList, toList)
import Data.Time (diffUTCTime, getCurrentTime)
import Data.Typeable (Typeable)
import Debug.Trace (traceShowId)
import GHC.IO.Unsafe (unsafePerformIO)
import GHC.TypeLits (KnownNat, Nat)
import HashedExpression.Internal.Expression
import HashedExpression.Internal.Normalize
import HashedExpression.Internal.Utils
import HashedExpression.Interp
import HashedExpression.Operation
import qualified HashedExpression.Operation
import HashedExpression.Prettify
import HashedExpression.Value
import Test.HUnit
import Test.Hspec
import Test.QuickCheck
import Var
import Prelude hiding ((^))

-- |
--
-- | Remove duplicate but also sort
removeDuplicate :: (Ord a) => [a] -> [a]
removeDuplicate = toList . fromList

-- |
vectorOfDifferent :: Eq a => Int -> Gen a -> Gen [a]
vectorOfDifferent sz gen = foldM f [] [1 .. sz]
  where
    f acc _ = (: acc) <$> gen `suchThat` (not . flip elem acc)

-- | Format
format :: [(String, String)] -> String
format = intercalate "\n" . map oneLine
  where
    oneLine (f, s) = f ++ ": " ++ s

-- |
inspect :: (Typeable d, Typeable rc) => Expression d rc -> Expression d rc
inspect x =
  unsafePerformIO $ do
    showExp x
    return x

-- | Vars list
type Vars = [[String]] -- Vars 0D, 1D, 2D, 3D, ..

mergeVars :: [Vars] -> Vars
mergeVars = foldl f [[], [], [], []]
  where
    f x y = map removeDuplicate $ zipWith (++) x y

genDouble :: Gen Double
genDouble = arbitrary `suchThat` inSmallRange
  where
    inSmallRange x = x >= 0 && x <= 10

-- |
genValMap ::
  forall size1D size2D1 size2D2.
  (KnownNat size1D, KnownNat size2D1, KnownNat size2D2) =>
  Vars ->
  Gen ValMaps
genValMap vars = do
  let sz1D = valueFromNat @size1D
      sz2D1 = valueFromNat @size2D1
      sz2D2 = valueFromNat @size2D2
  let [names0d, names1d, names2d, names3d] = vars
  list0d <- vectorOf (length names0d) genDouble
  let vm0 = Map.fromList . zip names0d $ map VScalar list0d
  list1d <- vectorOf (length names1d) . vectorOf sz1D $ genDouble
  let vm1 =
        Map.fromList . zip names1d . map (V1D . listArray (0, sz1D - 1)) $
          list1d
  list2d <- vectorOf (length names2d) . vectorOf (sz2D1 * sz2D2) $ genDouble
  let vm2 =
        Map.fromList
          . zip names2d
          . map (V2D . listArray ((0, 0), (sz2D1 - 1, sz2D2 - 1)))
          $ list2d
  let vm3 = Map.empty -- TODO: not testing with 3D yet.
  return $ Map.unions [vm0, vm1, vm2, vm3]

shouldApprox :: (HasCallStack, Approximable a) => a -> a -> Expectation
shouldApprox x y = assertBool msg (x ~= y)
  where
    msg = "Expected: " ++ prettifyShow y ++ "\nGot: " ++ prettifyShow x

infix 1 `shouldApprox`

liftE1 ::
  (Expression d1 et1 -> Expression d2 et2) ->
  (Expression d1 et1, Vars) ->
  (Expression d2 et2, Vars)
liftE1 op (e, v) = (op e, v)

liftE2 ::
  (Expression d1 et1 -> Expression d2 et2 -> Expression d3 et3) ->
  (Expression d1 et1, Vars) ->
  (Expression d2 et2, Vars) ->
  (Expression d3 et3, Vars)
liftE2 op (e1, v1) (e2, v2) = ((op e1 e2), (mergeVars [v1, v2]))

-------------------------------------------------------------------------------
primitiveScalarR :: Gen (Expression Scalar R, Vars)
primitiveScalarR = do
  name <- elements . map pure $ ['a' .. 'z']
  dbl <- genDouble
  r <- (`mod` 10) <$> arbitrary
  if r < 6
    then return ((variable name), [[(name)], [], [], []])
    else return ((constant dbl), [[], [], [], []])

primitiveScalarC :: Gen (Expression Scalar C, Vars)
primitiveScalarC = liftE2 (+:) <$> primitiveScalarR <*> primitiveScalarR

-------------------------------------------------------------------------------
primitive1DR ::
  forall n.
  KnownNat n =>
  Gen (Expression n R, Vars)
primitive1DR = do
  name <- elements . map ((++ "1") . pure) $ ['a' .. 'z']
  dbl <- genDouble
  r <- (`mod` 10) <$> arbitrary
  if r < 6
    then return ((variable1D @n name), [[], [(name)], [], []])
    else return ((constant1D @n dbl), [[], [], [], []])

primitive1DC ::
  forall n.
  KnownNat n =>
  Gen (Expression n C, Vars)
primitive1DC = liftE2 (+:) <$> primitive1DR @n <*> primitive1DR @n

-------------------------------------------------------------------------------
primitive2DR ::
  forall m n.
  (KnownNat m, KnownNat n) =>
  Gen (Expression '(m, n) R, Vars)
primitive2DR = do
  name <- elements . map ((++ "2") . pure) $ ['a' .. 'z']
  dbl <- genDouble
  r <- (`mod` 10) <$> arbitrary
  if r < 6
    then return (variable2D @m @n name, [[], [], [(name)], []])
    else return (constant2D @m @n dbl, [[], [], [], []])

primitive2DC ::
  forall m n.
  (KnownNat m, KnownNat n) =>
  Gen (Expression '(m, n) C, Vars)
primitive2DC = liftE2 (+:) <$> primitive2DR @m @n <*> primitive2DR @m @n

--type  = 10
-------------------------------------------------------------------------------
genScalarR ::
  forall default1D default2D1 default2D2.
  (KnownNat default1D, KnownNat default2D1, KnownNat default2D2) =>
  Int ->
  Gen (Expression Scalar R, Vars)
genScalarR size
  | size == 0 = primitiveScalarR
  | otherwise =
    let sub = genScalarR @default1D @default2D1 @default2D2 (size `div` 10)
        subC = genScalarC @default1D @default2D1 @default2D2 (size `div` 10)
        sub1D = gen1DR @default1D @default2D1 @default2D2 (size `div` 10)
        sub2D = gen2DR @default1D @default2D1 @default2D2 (size `div` 10)
        fromPiecewise = do
          numBranches <- elements [2 .. 4]
          branches <- vectorOf numBranches sub
          condition <- sub
          marks <- sort <$> vectorOfDifferent (numBranches - 1) arbitrary
          let vars = mergeVars $ map snd branches ++ [snd condition]
              exp = piecewise marks (fst condition) $ map fst branches
          return (exp, vars)
        binary op = liftE2 op <$> sub <*> sub
        unary op = liftE1 op <$> sub
     in oneof
          [ fromPiecewise,
            binary (+),
            binary (*),
            binary (*.),
            binary (-),
            binary (<.>),
            unary negate,
            unary (^ 2),
            liftE1 xRe <$> subC,
            liftE1 xIm <$> subC,
            liftE2 (<.>) <$> sub1D <*> sub1D,
            liftE2 (<.>) <$> sub2D <*> sub2D
          ]

-------------------------------------------------------------------------------
genScalarC ::
  forall default1D default2D1 default2D2.
  (KnownNat default1D, KnownNat default2D1, KnownNat default2D2) =>
  Int ->
  Gen (Expression Scalar C, Vars)
genScalarC size
  | size == 0 = primitiveScalarC
  | otherwise =
    let sub = genScalarC @default1D @default2D1 @default2D2 (size `div` 10)
        subR = genScalarR @default1D @default2D1 @default2D2 (size `div` 10)
        sub1D = gen1DC @default1D @default2D1 @default2D2 (size `div` 10)
        sub2D = gen2DC @default1D @default2D1 @default2D2 (size `div` 10)
        fromPiecewise = do
          numBranches <- elements [2 .. 4]
          branches <- vectorOf numBranches sub
          condition <- subR
          marks <- sort <$> vectorOfDifferent (numBranches - 1) arbitrary
          let vars = mergeVars $ map snd branches ++ [snd condition]
              exp = piecewise marks (fst condition) $ map fst branches
          return (exp, vars)
        binary op = liftE2 op <$> sub <*> sub
        unary op = liftE1 op <$> sub
     in oneof
          [ fromPiecewise,
            binary (+),
            binary (*),
            binary (*.),
            binary (-),
            binary (<.>),
            unary negate,
            unary (^ 2),
            liftE2 (+:) <$> subR <*> subR,
            liftE2 (<.>) <$> sub1D <*> sub1D,
            liftE2 (<.>) <$> sub2D <*> sub2D
          ]

-------------------------------------------------------------------------------
gen1DR ::
  forall n default2D1 default2D2.
  (KnownNat n, KnownNat default2D1, KnownNat default2D2) =>
  Int ->
  Gen (Expression n R, Vars)
gen1DR size
  | size == 0 = primitive1DR
  | otherwise =
    let sub = gen1DR @n @default2D1 @default2D2 (size `div` 10)
        subC = gen1DC @n @default2D1 @default2D2 (size `div` 10)
        subScalar = genScalarR @n @default2D1 @default2D2 (size `div` 10)
        fromPiecewise = do
          numBranches <- elements [2 .. 4]
          branches <- vectorOf numBranches sub
          condition <- sub
          marks <- sort <$> vectorOfDifferent (numBranches - 1) arbitrary
          let vars = mergeVars $ map snd branches ++ [snd condition]
              exp = piecewise marks (fst condition) $ map fst branches
          return (exp, vars)
        fromRotate = do
          amount <- elements [- (valueFromNat @n) .. valueFromNat @n]
          liftE1 (rotate amount) <$> sub
        binary op = liftE2 op <$> sub <*> sub
        unary op = liftE1 op <$> sub
     in oneof
          [ fromPiecewise,
            binary (+),
            binary (*),
            binary (-),
            unary negate,
            unary (^ 2),
            liftE1 xRe <$> subC,
            liftE1 xIm <$> subC,
            liftE2 (*.) <$> subScalar <*> sub,
            fromRotate,
            liftE1 (xRe . ft) <$> sub,
            liftE1 (xIm . ft) <$> sub
          ]

-------------------------------------------------------------------------------
gen1DC ::
  forall n default2D1 default2D2.
  (KnownNat n, KnownNat default2D1, KnownNat default2D2) =>
  Int ->
  Gen (Expression n C, Vars)
gen1DC size
  | size == 0 = primitive1DC
  | otherwise =
    let sub = gen1DC @n @default2D1 @default2D2 (size `div` 10)
        subR = gen1DR @n @default2D1 @default2D2 (size `div` 10)
        subScalar = genScalarC @n @default2D1 @default2D2 (size `div` 10)
        fromPiecewise = do
          numBranches <- elements [2 .. 4]
          branches <- vectorOf numBranches sub
          condition <- subR
          marks <- sort <$> vectorOfDifferent (numBranches - 1) arbitrary
          let vars = mergeVars $ map snd branches ++ [snd condition]
              exp = piecewise marks (fst condition) $ map fst branches
          return (exp, vars)
        fromRotate = do
          amount <- elements [- (valueFromNat @n) .. valueFromNat @n]
          liftE1 (rotate amount) <$> sub
        binary op = liftE2 op <$> sub <*> sub
        unary op = liftE1 op <$> sub
     in oneof
          [ fromPiecewise,
            binary (+),
            binary (*),
            binary (-),
            unary negate,
            unary (^ 2),
            liftE2 (+:) <$> subR <*> subR,
            liftE2 (*.) <$> subScalar <*> sub,
            liftE1 ft <$> sub,
            fromRotate
          ]

-------------------------------------------------------------------------------
gen2DR ::
  forall default1D m n.
  (KnownNat default1D, KnownNat m, KnownNat n) =>
  Int ->
  Gen (Expression '(m, n) R, Vars)
gen2DR size
  | size == 0 = primitive2DR
  | otherwise =
    let sub = gen2DR @default1D @m @n (size `div` 10)
        subC = gen2DC @default1D @m @n (size `div` 10)
        subScalar = genScalarR @default1D @m @n (size `div` 10)
        fromPiecewise = do
          numBranches <- elements [2 .. 4]
          branches <- vectorOf numBranches sub
          condition <- sub
          marks <- sort <$> vectorOfDifferent (numBranches - 1) arbitrary
          let vars = mergeVars $ map snd branches ++ [snd condition]
              exp = piecewise marks (fst condition) $ map fst branches
          return (exp, vars)
        fromRotate = do
          amount1 <- elements [- (valueFromNat @m) .. valueFromNat @m]
          amount2 <- elements [- (valueFromNat @n) .. valueFromNat @n]
          liftE1 (rotate (amount1, amount2)) <$> sub
        binary op = liftE2 op <$> sub <*> sub
        unary op = liftE1 op <$> sub
     in oneof
          [ fromPiecewise,
            binary (+),
            binary (*),
            binary (-),
            unary negate,
            unary (^ 2),
            liftE1 xRe <$> subC,
            liftE1 xIm <$> subC,
            liftE2 (*.) <$> subScalar <*> sub,
            fromRotate,
            liftE1 (xRe . ft) <$> sub,
            liftE1 (xIm . ft) <$> sub
          ]

-------------------------------------------------------------------------------
gen2DC ::
  forall default1D m n.
  (KnownNat default1D, KnownNat m, KnownNat n) =>
  Int ->
  Gen (Expression '(m, n) C, Vars)
gen2DC size
  | size == 0 = primitive2DC
  | otherwise =
    let sub = gen2DC @default1D @m @n (size `div` 10)
        subR = gen2DR @default1D @m @n (size `div` 10)
        subScalar = genScalarC @default1D @m @n (size `div` 10)
        fromPiecewise = do
          numBranches <- elements [2 .. 4]
          branches <- vectorOf numBranches sub
          condition <- subR
          marks <- sort <$> vectorOfDifferent (numBranches - 1) arbitrary
          let vars = mergeVars $ map snd branches ++ [snd condition]
              exp = piecewise marks (fst condition) $ map fst branches
          return (exp, vars)
        fromRotate = do
          amount1 <- elements [- (valueFromNat @m) .. valueFromNat @m]
          amount2 <- elements [- (valueFromNat @n) .. valueFromNat @n]
          liftE1 (rotate (amount1, amount2)) <$> sub
        binary op = liftE2 op <$> sub <*> sub
        unary op = liftE1 op <$> sub
     in oneof
          [ fromPiecewise,
            binary (+),
            binary (*),
            binary (-),
            unary negate,
            unary (^ 2),
            liftE2 (+:) <$> subR <*> subR,
            liftE2 (*.) <$> subScalar <*> sub,
            liftE1 ft <$> sub,
            fromRotate
          ]

-------------------------------------------------------------------------------
data Suite (size1D :: Nat) (size2D1 :: Nat) (size2D2 :: Nat) d et
  = Suite (Expression d et) ValMaps
  deriving (Show)

-------------------------------------------------------------------------------
type TestSuite = Suite Default1D Default2D1 Default2D2

type SuiteScalarR = TestSuite Scalar R

type SuiteScalarC = TestSuite Scalar C

type SuiteOneR = TestSuite Default1D R

type SuiteOneC = TestSuite Default1D C

type SuiteTwoR = TestSuite '(Default2D1, Default2D2) R

type SuiteTwoC = TestSuite '(Default2D1, Default2D2) C

-------------------------------------------------------------------------------
instance
  (KnownNat size1D, KnownNat size2D1, KnownNat size2D2) =>
  Arbitrary (Suite size1D size2D1 size2D2 Scalar R)
  where
  arbitrary = do
    (exp, vars) <- sized $ genScalarR @size1D @size2D1 @size2D2
    valMaps <- genValMap @size1D @size2D1 @size2D2 vars
    return $ Suite exp valMaps

instance
  (KnownNat size1D, KnownNat size2D1, KnownNat size2D2) =>
  Arbitrary (Suite size1D size2D1 size2D2 Scalar C)
  where
  arbitrary = do
    (exp, vars) <- sized $ genScalarC @size1D @size2D1 @size2D2
    valMaps <- genValMap @size1D @size2D1 @size2D2 vars
    return $ Suite exp valMaps

instance
  (KnownNat size1D, KnownNat size2D1, KnownNat size2D2) =>
  Arbitrary (Suite size1D size2D1 size2D2 size1D R)
  where
  arbitrary = do
    (exp, vars) <- sized $ gen1DR @size1D @size2D1 @size2D2
    valMaps <- genValMap @size1D @size2D1 @size2D2 vars
    return $ Suite exp valMaps

instance
  (KnownNat size1D, KnownNat size2D1, KnownNat size2D2) =>
  Arbitrary (Suite size1D size2D1 size2D2 size1D C)
  where
  arbitrary = do
    (exp, vars) <- sized $ gen1DC @size1D @size2D1 @size2D2
    valMaps <- genValMap @size1D @size2D1 @size2D2 vars
    return $ Suite exp valMaps

instance
  (KnownNat size1D, KnownNat size2D1, KnownNat size2D2) =>
  Arbitrary (Suite size1D size2D1 size2D2 '(size2D1, size2D2) R)
  where
  arbitrary = do
    (exp, vars) <- sized $ gen2DR @size1D @size2D1 @size2D2
    valMaps <- genValMap @size1D @size2D1 @size2D2 vars
    return $ Suite exp valMaps

instance
  (KnownNat size1D, KnownNat size2D1, KnownNat size2D2) =>
  Arbitrary (Suite size1D size2D1 size2D2 '(size2D1, size2D2) C)
  where
  arbitrary = do
    (exp, vars) <- sized $ gen2DC @size1D @size2D1 @size2D2
    valMaps <- genValMap @size1D @size2D1 @size2D2 vars
    return $ Suite exp valMaps

-------------------------------------------------------------------------------
instance Arbitrary (Expression Scalar R) where
  arbitrary = fst <$> sized (genScalarR @Default1D @Default2D1 @Default2D2)

instance Arbitrary (Expression Scalar C) where
  arbitrary = fst <$> sized (genScalarC @Default1D @Default2D1 @Default2D2)

instance KnownNat n => Arbitrary (Expression n R) where
  arbitrary = fst <$> sized (gen1DR @n @Default2D1 @Default2D2)

instance KnownNat n => Arbitrary (Expression n C) where
  arbitrary = fst <$> sized (gen1DC @n @Default2D1 @Default2D2)

instance (KnownNat m, KnownNat n) => Arbitrary (Expression '(m, n) R) where
  arbitrary = fst <$> sized (gen2DR @Default1D @m @n)

instance (KnownNat m, KnownNat n) => Arbitrary (Expression '(m, n) C) where
  arbitrary = fst <$> sized (gen2DC @Default1D @m @n)

-------------------------------------------------------------------------------
data ArbitraryExpresion = forall d et. (DimensionType d, ElementType et) => ArbitraryExpresion (Expression d et)

instance Show ArbitraryExpresion where
  show (ArbitraryExpresion exp) = show exp

instance Arbitrary ArbitraryExpresion where
  arbitrary =
    let option1 =
          fmap ArbitraryExpresion (arbitrary :: Gen (Expression Scalar R))
        option2 =
          fmap ArbitraryExpresion (arbitrary :: Gen (Expression Scalar C))
        option3 =
          fmap
            ArbitraryExpresion
            (arbitrary :: Gen (Expression Default1D R))
        option4 =
          fmap
            ArbitraryExpresion
            (arbitrary :: Gen (Expression Default1D C))
        option5 =
          fmap
            ArbitraryExpresion
            (arbitrary :: Gen (Expression '(Default2D1, Default2D2) R))
        option6 =
          fmap
            ArbitraryExpresion
            (arbitrary :: Gen (Expression '(Default2D1, Default2D2) C))
     in oneof [option1, option2, option3, option4, option5, option6]

-- |
getWrappedExp :: ArbitraryExpresion -> (ExpressionMap, NodeID)
getWrappedExp (ArbitraryExpresion (Expression n mp)) = (mp, n)

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

-- |
sz :: Expression d et -> Int
sz = IM.size . exMap

infix 1 `shouldNormalizeTo`

shouldNormalizeTo ::
  (HasCallStack, DimensionType d, ElementType et, Typeable et, Typeable d) =>
  Expression d et ->
  Expression d et ->
  IO ()
shouldNormalizeTo exp1 exp2 = do
  prettify (normalize exp1) `shouldBe` prettify (normalize exp2)
  normalize exp1 `shouldBe` normalize exp2

max_stars_repo_name: Turboscient/HashedExpression

max_issues_repo_name: Turboscient/HashedExpression

max_forks_repo_name: Turboscient/HashedExpression