text:
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE FlexibleContexts #-}

module Main where
import Control.Monad.ST
import Control.Monad.Trans.Reader
import System.Random.MWC
import Data.Word
import Data.Maybe (isJust)
import Control.Monad (replicateM)
import Control.Monad.State
import System.Exit
import Data.List (sort)
import qualified Data.Vector.Generic hiding (replicateM, sum, product)
import Statistics.Distribution
import Statistics.Distribution.Exponential
import Graphics.Rendering.Chart
import Graphics.Rendering.Chart.Backend.Diagrams
import Data.Default.Class
import Control.Lens

-- Ereignis Objekt
data Event = Arrival {
    time :: Double,
    customersWaiting :: Int
  } | Departure {
    time :: Double,
    customersWaiting :: Int,
    waitingTime :: Double
  }deriving (Show)

data Simulation = Simulation [Double] Double (Maybe Double)

-- exponentiell verteilten Zufallswert für die Bedienzeit
rndServiceTime :: Rand Double
rndServiceTime = genContV $ exponential 0.9

-- Zeit bis zur nächsten Ankunft eines Kunden
rndTimeNextCustomer :: Rand Double
rndTimeNextCustomer = genContV $ exponential 0.6 -- (1/0.6)

departure :: Int -> Simulation -> IO [Event]
departure n (Simulation [] _ _) = exitWith (ExitFailure 1)
departure n (Simulation _ _ Nothing) = exitWith (ExitFailure 1)
departure n (Simulation (wc:ws) nc (Just t)) = do
    let l = length ws -- get amount of waiting customers
    st <- if l == 0 then return Nothing else Just <$> (+) t <$> runRandIO rndServiceTime -- service time of the arrived customer
    (Departure t l (t-wc):) <$> step (n-1) (Simulation ws nc st) -- create the event

arrival :: Int -> Simulation -> IO [Event]
arrival n (Simulation wc t st) = do
    nc <- (+) t <$> runRandIO rndTimeNextCustomer  -- when the next customer will arrive
    let dp = wc ++ [t] -- times of waiting customers
    st <- if isJust st then return st else Just <$> (+) t <$> runRandIO rndServiceTime  -- service time of the arrived customer if there is not already a customer served
    (Arrival t (length dp):) <$> step (n-1) (Simulation dp nc st) -- create the event

isArrival :: Simulation -> Bool
isArrival (Simulation [] _ _) = True
isArrival (Simulation _ _ Nothing) = True
isArrival (Simulation _ ac (Just dc)) = ac <= dc

step :: Int -> Simulation -> IO [Event]
step 0 _ = return []
step n sim = if isArrival sim
             then arrival n sim
             else departure n sim

runSim :: Int -> IO [Event]
runSim n = step n (Simulation [] 0 Nothing)

main = do
    events <-  runSim 100
    let xy = map (\e -> (time e, customersWaiting e)) events
    let formattedData = concat $ map format $ zipWith (\a b -> (a, fst b, snd b)) ((head xy):(init xy)) $ zip xy ((tail xy) ++ [last xy])
    renderableToFile def "simulation.svg" $ chart $ formattedData
    putStrLn $ ("Maximale Wartezeit: " ++) . show $ maximum $ [ w | x@(Departure _ _ w) <- events ]
    putStrLn $ ("Maximale Anzahl wartender Kunden: " ++) . show $ maximum $ map customersWaiting events
  where
    format ((t1, cw1), (t2, cw2), (t3, cw3)) = [(average t1 t2, cw1), (average t1 t2, cw2), (average t2 t3, cw2)]
    average x y = (x+y)/2

type Rand0 s a = ReaderT (Gen s) (ST s) a
type Rand a = (forall s. Rand0 s a)  -- the random-sampling monad

-- Draw from the continuous distribution d
genContV d = ask >>= genContVar d

-- Provide a seed for the PRNG and return a draw from the random sampler
runRandV :: Data.Vector.Generic.Vector v Word32 => Rand a -> v Word32 -> a
runRandV rand seeds =
  runST $ initialize seeds >>= runReaderT rand

-- Seed the PRNG with data from the system's fast source of pseudo-random numbers,
-- then return a draw from the random sampler
runRandIO :: Rand a -> IO a
runRandIO rand = do
  gen <- createSystemRandom
  seeds <- fromSeed <$> save gen
  return $ runRandV rand seeds

-- Render plot

chart xs = toRenderable layout
  where
    plot = plot_lines_values .~ [xs]
             $ def
    layout = layout_title .~ "Wartende Kunden"
           $ layout_x_axis . laxis_title .~ "Zeit"
           $ layout_y_axis . laxis_title .~ "Kunden"
           $ layout_plots .~ [toPlot plot]
           $ def

max_stars_repo_name: Zortaniac/simulation

max_issues_repo_name: Zortaniac/simulation

max_forks_repo_name: Zortaniac/simulation


text:
module Reexports
  ( Matrix
  , Vector
  , R
  , toLists
  , toList
  , matrix
  , fromLists
  , fromList
  ) where

import Numeric.LinearAlgebra.Data
  ( Matrix
  , R
  , Vector
  , fromList
  , fromLists
  , matrix
  , toList
  , toLists
  )

ext: hs

max_stars_repo_name: sebashack/mlToolBox

max_issues_repo_name: sebashack/mlToolBox

max_forks_repo_name: sebashack/mlToolBox


text:
------------------------------------------------
-- |
-- Module    :  Numeric.Matrix.Integral.SmithNF
-- Copyright :  (c) Jun Yoshida 2019
-- License   :  BSD3
--
-- Compute Smith normal form
--
------------------------------------------------

module Numeric.Matrix.Integral.SmithNF where

import Control.Monad (when, unless, forM_, guard, join)
import Control.Monad.ST
import Control.Monad.Loops (whileM_, whileJust_)

import Data.List as L
import Data.STRef

import Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra.Devel

--import Numeric.Matrix.Integral.HNFLLL (hnfLLLST)
import Numeric.Matrix.Integral.NormalForms (hermiteNFST)

-- | Find off-diagonal entries and calculate min of min{row,col} among them
findDiagSz :: LA.Matrix LA.Z -> Maybe Int
findDiagSz ma =
  let offdiag = filter (uncurry (/=)) $ LA.find (/=0) ma
  in L.foldl' bin Nothing offdiag
  where
    bin Nothing (i,j) = Just (min i j)
    bin (Just k) (i,j) = Just $ L.foldl' min k [i,j]

-- | Compute "pre-Smith" form with a function calculating Hermite normal forms
-- "Pre-Smith" means diag(c_1,c_2,..) with not necessarily c_i|c_{i+1}
preSmithNFST :: STMatrix s LA.Z -> STMatrix s LA.Z ->  STMatrix s LA.Z -> ST s ()
preSmithNFST stMatUL stMatA stMatURt = do
  whileJust_ (liftSTMatrix findDiagSz stMatA) $ \diagSz -> do
    -- Extract a submatrix containing all the off-diagonal entries.
    stMatAex <- thawMatrix =<< extractMatrix stMatA (FromRow diagSz) (FromCol diagSz)
    stMatULex <- thawMatrix =<< extractMatrix stMatUL (FromRow diagSz) AllCols
    stMatURtex <- thawMatrix =<< extractMatrix stMatURt (FromRow diagSz) AllCols
    -- Compute the Hermite normal forms of the extracted submatrix and of the transpose of the result.
    --hnfLLLST stMatULex stMatAex
    join $ hermiteNFST <$> newMatrix 0 0 0 <*> pure stMatULex <*> pure stMatAex
    stMatAext <- thawMatrix =<< (LA.tr' <$> freezeMatrix stMatAex)
    --hnfLLLST stMatURtex stMatAext
    join $ hermiteNFST <$> newMatrix 0 0 0 <*> pure stMatURtex <*> pure stMatAext
    -- Write out the result
    setMatrix stMatA diagSz diagSz =<< (LA.tr' <$> freezeMatrix stMatAext)
    setMatrix stMatUL diagSz 0 =<< freezeMatrix stMatULex
    setMatrix stMatURt diagSz 0 =<< freezeMatrix stMatURtex

isHeadGCD :: LA.Vector LA.Z -> Bool
isHeadGCD = isHeadGCD' . LA.toList
  where
    isHeadGCD' [] = True
    isHeadGCD' (x:xs) = L.foldl' (\b y -> b && rem y x == 0) True xs

-- Compute the submatrix containing all the non-zero diagonal entries
extractMaxNZDiag :: LA.Matrix LA.Z -> (Int,LA.Matrix LA.Z)
extractMaxNZDiag mx = runST $ do
  dRef <- newSTRef 0
  let p d = d < uncurry min (LA.size mx) && (mx!d!d) /= 0
  whileM_ (p <$> readSTRef dRef) $ modifySTRef' dRef (+1)
  d <- readSTRef dRef
  return (d,LA.subMatrix (0,0) (d,d) mx)

-- | Normalize diagonals to form a divisible chain; i.e. c_1 | c_2 | c_3 | ...
normalizeST :: STMatrix s LA.Z -> STMatrix s LA.Z ->  STMatrix s LA.Z -> ST s ()
normalizeST stMatUL stMatA stMatURt = do
  (dlen,matD) <- liftSTMatrix extractMaxNZDiag stMatA
  stMatD <- thawMatrix matD
  stMatULex <- thawMatrix =<< extractMatrix stMatUL (RowRange 0 (dlen-1)) AllCols
  stMatURtex <- thawMatrix =<< extractMatrix stMatURt (RowRange 0 (dlen-1)) AllCols
  normalizeST' stMatULex stMatD stMatURtex
  setMatrix stMatA 0 0 =<< freezeMatrix stMatD
  setMatrix stMatUL 0 0 =<< freezeMatrix stMatULex
  setMatrix stMatURt 0 0 =<< freezeMatrix stMatURtex
  where
    normalizeST' stMatULex stMatD stMatURtex = do
      b <- (isHeadGCD . LA.takeDiag) <$> freezeMatrix stMatD
      n <- liftSTMatrix LA.rows stMatD
      unless b $ do
        forM_ [1..(n-1)] $ \i -> do
          readMatrix stMatD i i >>= writeMatrix stMatD i 0
          rowOper (AXPY 1 i 0 AllCols) stMatURtex
        preSmithNFST stMatULex stMatD stMatURtex
      when (n > 1) $ do
        stMatDex <- thawMatrix =<< extractMatrix stMatD (FromRow 1) (FromCol 1)
        stMatULexex <- thawMatrix =<< extractMatrix stMatULex (FromRow 1) AllCols
        stMatURtexex <- thawMatrix =<< extractMatrix stMatURtex (FromRow 1) AllCols
        normalizeST' stMatULexex stMatDex stMatURtexex
        setMatrix stMatD 1 1 =<< freezeMatrix stMatDex
        setMatrix stMatULex 1 0 =<< freezeMatrix stMatULexex
        setMatrix stMatURtex 1 0 =<< freezeMatrix stMatURtexex

-- | Compute Smith normal form together with transform unimodular matrices.
smithNF :: LA.Matrix LA.Z -> (LA.Matrix LA.Z, LA.Matrix LA.Z, LA.Matrix LA.Z)
smithNF matA = runST $ do
  stMatA <- thawMatrix matA
  stMatUL <- thawMatrix $ LA.ident (LA.rows matA)
  stMatURt <- thawMatrix $ LA.ident (LA.cols matA)
  preSmithNFST stMatUL stMatA stMatURt
  normalizeST stMatUL stMatA stMatURt
  resMatS <- freezeMatrix stMatA
  resMatUL <- freezeMatrix stMatUL
  resMatURt <- freezeMatrix stMatURt
  return (resMatUL, resMatS, LA.tr' resMatURt)

-- | For debug.
-- Compute "pre-Smith form."
preSmithNF :: LA.Matrix LA.Z -> (LA.Matrix LA.Z, LA.Matrix LA.Z, LA.Matrix LA.Z)
preSmithNF matA = runST $ do
  stMatA <- thawMatrix matA
  stMatUL <- thawMatrix $ LA.ident (LA.rows matA)
  stMatURt <- thawMatrix $ LA.ident (LA.cols matA)
  preSmithNFST stMatUL stMatA stMatURt
  resMatS <- freezeMatrix stMatA
  resMatUL <- freezeMatrix stMatUL
  resMatURt <- freezeMatrix stMatURt
  return (resMatUL, resMatS, LA.tr' resMatURt)

max_stars_repo_name: Junology/linkedgram

max_issues_repo_name: Junology/linkedgram

max_forks_repo_name: Junology/linkedgram


text:
module NaiveBayes.GaussianNB
  ( train
  , predict
  , GaussianNB(..)
  ) where
import           Data.List           ((++))
import qualified Data.List           as L
import           Data.Map            ((!))
import qualified Data.Map            as M
import qualified Data.Vector.Unboxed as U
import           Prelude             hiding ((++))
import           Statistics.Matrix   as S
import           Statistics.Sample   as S

-- Implementation taken from
-- https://chrisalbon.com/machine_learning/naive_bayes/naive_bayes_classifier_from_scratch/

-- | Usage Example
--
-- trainFeatures :: S.Matrix
-- trainFeatures = S.fromRowLists
--   [ [6.00, 180, 12]
--   , [5.92, 190, 11]
--   , [5.58, 170, 12]
--   , [5.92, 165, 10]
--   , [5.00, 100, 6]
--   , [5.50, 150, 8]
--   , [5.42, 130, 7]
--   , [5.75, 150, 9]
--   ]
--
-- trainLabels :: U.Vector Int
-- trainLabels = U.fromList [0, 0, 0, 0, 1, 1, 1, 1]
--
-- trainedData = train trainFeatures trainLabels
-- predict trainedData [6.00, 130, 8]
-- 1

data GaussianNB = GaussianNB
  { cMeans      :: M.Map Int [Double]
  , cVariances  :: M.Map Int [Double]
  , cPriorProbs :: M.Map Int Double
  }

  deriving (Show)

labelsCounts :: U.Vector Int -> M.Map Int Int
labelsCounts = U.foldr' (\x acc -> M.insertWith (+) x 1 acc) M.empty

labelsFeatures :: S.Matrix -> U.Vector Int -> M.Map Int [U.Vector Double]
labelsFeatures features labels =
  L.foldr
    (\(l, fs) acc -> M.insertWith (++) l [fs] acc)
    M.empty
    (zip' labels features)

zip' :: U.Vector Int -> S.Matrix -> [(Int, U.Vector Double)]
zip' vec mx = U.ifoldr' (\i x acc -> (x, S.row mx i):acc) [] vec

means :: M.Map Int [U.Vector Double] -> M.Map Int [Double]
means = fmap (fmap S.mean . S.toRows . S.transpose . S.fromRows)

variances :: M.Map Int [U.Vector Double] -> M.Map Int [Double]
variances =
  fmap (fmap S.variance . S.toRows . S.transpose . S.fromRows)

priorProbs :: M.Map Int Int -> Int -> M.Map Int Double
priorProbs mp totalCount =
  fmap (\x -> fromIntegral x / fromIntegral totalCount) mp

train :: Matrix -> U.Vector Int -> GaussianNB
train features labels =
  let countMap = labelsCounts labels
      featuresMap = labelsFeatures features labels
      meansMap = means featuresMap
      variancesMap = variances featuresMap
      priorProbsMap = priorProbs countMap (U.length labels)
  in  GaussianNB meansMap variancesMap priorProbsMap

predict :: GaussianNB -> [Double] -> Int
predict cf feature =
  fst $ L.maximumBy
    (\(_, a) (_, b) -> a `compare` b)
    (M.toList (posteriorProbs cf feature))

posteriorProbs :: GaussianNB -> [Double] -> M.Map Int Double
posteriorProbs (GaussianNB ms vs pp) feature =
  M.mapWithKey (\k x -> x * likelihoods k) pp
    where likelihoods k = L.product $
            L.map
              (\(x, m, v) -> exp ((- ((x - m) ** 2)) / (2 * v)) / (sqrt (2 * pi * v)))
              (L.zip3 feature (ms ! k) (vs ! k))

max_stars_repo_name: abs-zero/haskellml

max_issues_repo_name: abs-zero/haskellml

max_forks_repo_name: abs-zero/haskellml


text:
module ZeroStuff (zeroStuff) where

import IQ
import Data.Complex
import Control.Monad
import Pipes

zeroStuff :: (Monad m, Num a) => Int -> Pipe a a m ()
zeroStuff factor = forever $ await >>= yield >> replicateM_ (factor - 1) (yield 0)
{-# INLINE zeroStuff #-}

max_stars_repo_name: hexagonal-sun/ayeQ

max_issues_repo_name: hexagonal-sun/ayeQ

max_forks_repo_name: hexagonal-sun/ayeQ


text:
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE EmptyCase #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE LambdaCase #-}

module Models.Integrals.Approximate2  where

import Algebra.Classes
import qualified Algebra.Morphism.Affine as A
import Prelude hiding (Num(..), Fractional(..), (^), product, sum, pi, sqrt
                      , exp)
import qualified Models.Field
import qualified Algebra.Linear.Chebyshev as Chebyshev
import Data.Complex
import qualified Algebra.Morphism.Polynomial.Multi as Multi
import Algebra.Morphism.Polynomial.Multi hiding (constPoly)
import qualified Algebra.Morphism.LinComb as LC
import Control.Applicative

import Models.Integrals.Types

type family Env γ where
  Env 'Unit = ()
  Env (a ':× b) = (Env a, Env b)
  Env 'R = Double

type RR = Double

newtype PointWise x a = PW (x -> a) deriving Functor
instance Applicative (PointWise x) where
  pure x = PW $  \_ -> x
  PW f <*> PW a = PW (\x -> (f x) (a x))

fromPointWise :: PointWise x a -> x -> a
fromPointWise (PW x) = x

instance Additive a => Additive (PointWise x a) where
  zero = pure zero
  (+) = liftA2 (+)

instance Group a => Group (PointWise x a) where
  negate = fmap negate
  (-) = liftA2 (-)

instance AbelianAdditive a => AbelianAdditive (PointWise x a)

instance Multiplicative a => Multiplicative (PointWise x a) where
  one  = pure one
  (*) = liftA2 (*)

instance Division a => Division (PointWise x a) where
  recip  = fmap recip
  (/) = liftA2 (/)

instance Roots a => Roots (PointWise x a) where
  root n  = fmap (root n)

instance Transcendental a => Transcendental (PointWise x a) where
  pi = pure pi
  -- log = fmap log
  -- sin = fmap sin
  -- cos = fmap cos
  -- asin = fmap asin
  -- acos = fmap acos
  -- atan = fmap atan
  -- sinh = fmap sinh
  -- cosh = fmap cosh
  -- asinh = fmap asinh
  -- acosh = fmap acosh
  -- atanh = fmap atanh
  exp = fmap exp

instance Multiplicative a => Scalable (PointWise x a) (PointWise x a) where
  (*^) = (*)

instance Ring a => Ring (PointWise x a) where
  fromInteger = PW . const . fromInteger

instance Field a => Field (PointWise x a) where
  fromRational = PW . const . fromRational

type S γ = PointWise (Env γ) RR

lk :: Var γ -> Env γ -> RR
lk Get (_,x) = x
lk (Weaken x) (ρ,_) = lk x ρ

lk' :: Var γ -> S γ
lk' v = PW $ lk v

evalC :: forall γ. Coef γ -> S γ
evalC (Coef c) = LC.eval (\x -> PW @(Env γ) $ \_ρ -> Models.Field.eval @RR x)
                         (exp . evalP) c

evalP :: forall γ. Poly γ -> S γ
evalP = Multi.eval evalC evalE

evalE :: forall γ. Elem γ -> S γ
evalE = \case
  Supremum dir es -> case dir of
    Min ->  foldr (liftA2 min) (pure ( 1/0)) (evalP <$> es)
    Max ->  foldr (liftA2 max) (pure (-1/0)) (evalP <$> es)
  Vari x -> lk' x
  CharFun x -> fmap (\y -> if y >= 0 then 1 else 0) (evalP x)

approxIntegrals :: Int -> Env γ -> P γ -> RR
approxIntegrals n ρ =
  let evP x = fromPointWise (evalP x) ρ
  in \case
      Add a b -> approxIntegrals n ρ a + approxIntegrals n ρ b
      Div a b -> approxIntegrals n ρ a / approxIntegrals n ρ b
      Integrate (mkSuprema -> (lo,hi)) e ->
        realPart (Chebyshev.integral @Double @C n (evP lo) (evP hi) $
                  \x -> approxIntegrals n (ρ,x) e :+ 0)
      Ret x -> evP x
      Cond (IsNegative c) e -> if A.eval Models.Field.eval (flip lk ρ) c <= 0 then approxIntegrals n ρ e else 0
      Cond (IsZero _) _ -> error "approxIntegrals: equality not eliminated?"


text:
-- | The method of moments can be used to estimate a number of commonly used distributions.  This module is still under construction as I work out the best way to handle morphisms from the Moments3 type to types of other distributions.  For more information, see the wikipedia entry: <https://en.wikipedia.org/wiki/Method_of_moments_(statistics)>

module HLearn.Models.Distributions.Univariate.Geometric
    ( Geometric
    )
    where

import Control.DeepSeq
import GHC.TypeLits
import qualified Data.Vector.Unboxed as U
import Data.Vector.Unboxed.Deriving

import qualified Statistics.Distribution as S
import Statistics.Distribution.Geometric

import HLearn.Algebra
import HLearn.Models.Distributions.Common
import HLearn.Models.Distributions.Univariate.Internal.Moments
import HLearn.Models.Distributions.Visualization.Gnuplot

-------------------------------------------------------------------------------
-- Geometric

newtype Geometric prob dp = Geometric {  moments :: (Moments3 dp) }
    deriving (Read,Show,Eq,Ord,Monoid,Group)

instance (Num dp) => HomTrainer (Geometric prob dp) where
    type Datapoint (Geometric prob dp) = dp
    train1dp dp = Geometric $ train1dp dp

instance (Num dp) => Probabilistic (Geometric prob dp) where
    type Probability (Geometric prob dp) = prob

instance (Integral dp, Floating prob) => PDF (Geometric prob dp) where
    pdf dist dp = p*(1-p)^^dp
        where
            p = geo_p dist

instance
    ( PDF (Geometric prob dp)
    , Show prob
    , Show dp
    , Ord dp
    , Ord prob
    , Fractional prob
    , RealFrac prob
    , Integral dp
    ) => PlottableDistribution (Geometric prob dp)
        where

    plotType _ = Points

    samplePoints dist = [min..max]
        where
            min = 0
            max = maximum [20,round $ 3*(fromIntegral $ mean dist)]

geo_p :: (Fractional prob, Integral dp) => Geometric prob dp -> prob
geo_p (Geometric dist) = 1/((fromIntegral $ m1 dist)/(fromIntegral $ m0 dist) +1)

instance (Integral dp, RealFrac prob) => Mean (Geometric prob dp) where
    mean dist = round $ 1/(geo_p dist)

instance (Integral dp, Fractional prob) => Variance (Geometric prob dp) where
    variance dist = (1-p)/(p*p)
        where
            p = geo_p dist

max_stars_repo_name: silky/HLearn

max_issues_repo_name: silky/HLearn

max_forks_repo_name: silky/HLearn


text:
{-# LANGUAGE DeriveFunctor, DeriveDataTypeable,DeriveGeneric  #-}
module Statistics.Test.Types (
    Test(..)
  , isSignificant
  , TestResult(..)
  , significant
  , PositionTest(..)
  ) where

import Control.DeepSeq  (NFData(..))
import Control.Monad    (liftM3)
import Data.Aeson       (FromJSON, ToJSON)
import Data.Binary      (Binary (..))
import Data.Data (Typeable, Data)
import GHC.Generics

import Statistics.Types (PValue)

-- | Result of hypothesis testing
data TestResult = Significant    -- ^ Null hypothesis should be rejected
                | NotSignificant -- ^ Data is compatible with hypothesis
                  deriving (Eq,Ord,Show,Typeable,Data,Generic)

instance Binary   TestResult where
  get = do
      sig <- get
      if sig then return Significant else return NotSignificant
  put = put . (== Significant)
instance FromJSON TestResult
instance ToJSON   TestResult
instance NFData   TestResult

-- | Result of statistical test.
data Test distr = Test
  { testSignificance :: !(PValue Double)
    -- ^ Probability of getting value of test statistics at least as
    --   extreme as measured.
  , testStatistics   :: !Double
    -- ^ Statistic used for test.
  , testDistribution :: distr
    -- ^ Distribution of test statistics if null hypothesis is correct.
  }
  deriving (Eq,Ord,Show,Typeable,Data,Generic,Functor)

instance (Binary   d) => Binary   (Test d) where
  get = liftM3 Test get get get
  put (Test sign stat distr) = put sign >> put stat >> put distr
instance (FromJSON d) => FromJSON (Test d)
instance (ToJSON   d) => ToJSON   (Test d)
instance (NFData   d) => NFData   (Test d) where
  rnf (Test _ _ a) = rnf a

-- | Check whether test is significant for given p-value.
isSignificant :: PValue Double -> Test d -> TestResult
isSignificant p t
  = significant $ p >= testSignificance t

-- | Test type for test which compare positional (mean,median etc.)
--   information of samples.
data PositionTest
  = SamplesDiffer
    -- ^ Test whether samples differ in position. Null hypothesis is
    --   samples are not different
  | AGreater
    -- ^ Test if first sample (A) is larger than second (B). Null
    --   hypothesis is first sample is not larger than second.
  | BGreater
    -- ^ Test if second sample is larger than first.
  deriving (Eq,Ord,Show,Typeable,Data,Generic)

instance Binary   PositionTest where
  get = do
    i <- get
    case (i :: Int) of
      0 -> return SamplesDiffer
      1 -> return AGreater
      2 -> return BGreater
      _ -> fail "Invalid PositionTest"
  put SamplesDiffer = put (0 :: Int)
  put AGreater      = put (1 :: Int)
  put BGreater      = put (2 :: Int)
instance FromJSON PositionTest
instance ToJSON   PositionTest
instance NFData   PositionTest

-- | significant if parameter is 'True', not significant otherwise
significant :: Bool -> TestResult
significant True  = Significant
significant False = NotSignificant


text:
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DeriveAnyClass #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE UnicodeSyntax #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE CPP #-}

module NQS.SR
  ( -- * Stochastic Reconfiguration
    --
    -- Globally, the SR algorithm looks more or less like this
    --
    -- @
    --     for n in {0,...,maxIter - 1} do
    --         (F, ∂ψ) <- sample ψ
    --         S <- covariance⁽ⁿ⁾ ∂ψ
    --         δ <- (S + λ⁽ⁿ⁾)⁻¹ F
    --         ψ <- ψ - learningRate * δ
    --     done
    -- @

    -- * AXPY
    -- * Sampling
    -- * Constructing S
    -- * Solving S
    sr
  , IterInfo(..)
  , HasIter(..)
  , HasSampler(..)
  , HasSolver(..)
  , HasMoments(..)
  , HasForceNorm(..)
  ) where

import           Prelude                 hiding ( zipWith
                                                , zipWithM
                                                , map
                                                , mapM
                                                )

import           GHC.Generics                   ( Generic )

import           Debug.Trace
import           Control.Exception              ( assert )
import           Control.Monad.Identity         ( Identity(..) )
import           Control.Monad.ST
import           System.IO.Unsafe               ( unsafePerformIO )
import           Foreign.Storable
import           Foreign.ForeignPtr
import           Data.Vector.Storable           ( Vector )
import qualified Data.Vector.Storable          as V
import           Data.Vector.Storable.Mutable   ( MVector )
import qualified Data.Vector.Storable.Mutable  as MV
import qualified Data.Vector.Unboxed

import           Data.Singletons
import           Data.Complex
import           Data.Semigroup                 ( (<>) )
import           Control.Monad                  ( (>=>) )
import           Control.Monad.Primitive

import           Control.DeepSeq
import           System.CPUTime
import Data.Aeson

import qualified NQS.CG                        as CG
import           NQS.CG                         ( Operator )
import           NQS.Rbm (Rbm(..))
import           NQS.Rbm.Mutable
import           NQS.Internal.BLAS
import           NQS.Internal.LAPACK
import           NQS.Internal.Types
import           NQS.Internal.Hamiltonian
import           NQS.Internal.Sampling
import           NQS.Internal.Rbm (unsafeFreezeRbm)

import           Lens.Micro
import           Lens.Micro.TH
import           Lens.Micro.Extras

import           GHC.Float                      ( int2Float )

import           Data.Aeson
import qualified Data.ByteString.Lazy.Char8     as BS

-- | In Stochastic Reconfiguration we only ever deal with 'Direct' vectors.
type V = MDenseVector 'Direct

-- | A shorter name for dense matrices.
type M orient s a = MDenseMatrix orient s a

data SolverStats = SolverStats
  { _solverStatsIters :: {-# UNPACK #-}!Int
  , _solverStatsErr   :: {-# UNPACK #-}!ℝ
  , _solverStatsTime  :: {-# UNPACK #-}!Double
  }

makeFields ''SolverStats

instance ToJSON SolverStats where
  toJSON stats =
    object [ "iters" .= (stats ^. iters)
           , "error" .= (stats ^. err)
           , "time"  .= (stats ^. time)
           ]
  toEncoding stats =
    pairs $ "iters" .= (stats ^. iters)
         <> "error" .= (stats ^. err)
         <> "time"  .= (stats ^. time)

instance FromJSON SolverStats where
  parseJSON = withObject "SolverStats" $ \v ->
    SolverStats <$> v .: "iters"
                <*> v .: "error"
                <*> v .: "time"

data SamplerStats = SamplerStats
  { _samplerStatsMoments  :: {-# UNPACK #-}!(Vector ℂ)
  , _samplerStatsStdDev   :: !(Maybe ℝ)
  , _samplerStatsDim      :: {-# UNPACK #-}!Int
  , _samplerStatsTime     :: {-# UNPACK #-}!Double
  }

makeFields ''SamplerStats

instance ToJSON SamplerStats where
  toJSON stats =
    object [ "moments" .= (stats ^. moments)
           , "stddev"  .= (stats ^. stdDev)
           , "dim"     .= (stats ^. dim)
           , "time"    .= (stats ^. time)
           ]
  toEncoding stats =
    pairs $ "moments" .= (stats ^. moments)
         <> "stddev"  .= (stats ^. stdDev)
         <> "dim"     .= (stats ^. dim)
         <> "time"    .= (stats ^. time)

instance FromJSON SamplerStats where
  parseJSON = withObject "SamplerStats" $ \v ->
    SamplerStats <$> v .: "moments"
                 <*> v .: "stddev"
                 <*> v .: "dim"
                 <*> v .: "time"

data IterInfo = IterInfo
  { _iterInfoIter      :: {-# UNPACK #-}!Int
  , _iterInfoState     :: Rbm
  , _iterInfoSampler   :: !SamplerStats
  , _iterInfoSolver    :: !SolverStats
  , _iterInfoForceNorm :: {-# UNPACK #-}!ℝ
  , _iterInfoDeltaNorm :: {-# UNPACK #-}!ℝ
  }

makeFields ''IterInfo

instance ToJSON IterInfo where
  toJSON stats =
    object [ "iter"      .= (stats ^. iter)
           , "state"     .= (stats ^. state)
           , "sampler"   .= (stats ^. sampler)
           , "solver"    .= (stats ^. solver)
           , "forceNorm" .= (stats ^. forceNorm)
           , "deltaNorm" .= (stats ^. deltaNorm)
           ]
  toEncoding stats =
    pairs $ "iter"      .= (stats ^. iter)
         <> "state"     .= (stats ^. state)
         <> "sampler"   .= (stats ^. sampler)
         <> "solver"    .= (stats ^. solver)
         <> "forceNorm" .= (stats ^. forceNorm)
         <> "deltaNorm" .= (stats ^. deltaNorm)

instance FromJSON IterInfo where
  parseJSON = withObject "IterInfo" $ \v ->
    IterInfo <$> v .: "iter"
             <*> v .: "state"
             <*> v .: "sampler"
             <*> v .: "solver"
             <*> v .: "forceNorm"
             <*> v .: "deltaNorm"

-- AXPY
-- -----------------------------------------------------------------------------

type AxpyFun v m = ℂ -> V (PrimState m) ℂ -> v (PrimState m) -> m ()

-- | \(Y \leftarrow \alpha X + Y\) for 'Column' matrices.
axpyMatrix
  :: (PrimMonad m, Storable a, AXPY (MDenseVector 'Direct) a)
  => a
  -> MDenseVector 'Direct (PrimState m) a
  -> MDenseMatrix 'Column (PrimState m) a
  -> m ()
axpyMatrix !α !x !y@(MDenseMatrix !rows !cols !stride !buff) =
  assert (rows * cols == x ^. dim) $ loop 0
 where
  getX !i = slice (i * rows) rows x
  getY !i = unsafeColumn i y
  loop !i | i < cols  = axpy α (getX i) (getY i) >> loop (i + 1)
          | otherwise = return ()

-- | Performs @y <- αx + y@ where @x@ is a vector and @y@ -- an RBM.
axpyRbm :: AxpyFun MRbm IO
axpyRbm α x y =
  let n = sizeVisible y
      m = sizeHidden y
  in  assert (x ^. dim == size y) $ do
        withVisible y $ \a -> do
          axpy α (slice 0 n x) a
          -- unsafeFreeze a >>= \a' -> print (a' ^. buffer)
        withHidden y $ \b -> do
          axpy α (slice n m x) b
          -- unsafeFreeze b >>= \b' -> print (b' ^. buffer)
        withWeights y $ \w -> do
          axpyMatrix α (slice (n + m) (n * m) x) w
          -- unsafeFreeze w >>= \w' -> print (w' ^. buffer)

-- Monte-Carlo sampling
-- -----------------------------------------------------------------------------

-- | Sampling
type SampleFun v m
  = v (PrimState m)
 -> V (PrimState m) ℂ
 -> M 'Row (PrimState m) ℂ
 -> m SamplerStats

-- | Monte-Carlo sampling an RBM
sampleRbm :: MCConfig -> Hamiltonian -> SampleFun MRbm IO
sampleRbm config hamiltonian = doSample
 where
  doSample rbm force derivatives = do
    t1                    <- getCPUTime
    moments               <- MV.new 2
    (dimension, variance) <- sampleGradients config
                                             hamiltonian
                                             rbm
                                             moments
                                             force
                                             derivatives
    t2 <- getCPUTime
    V.unsafeFreeze moments >>= \moments' -> return
      (SamplerStats moments'
                    (sqrt <$> variance)
                    dimension
                    (fromIntegral (t2 - t1) * 1.0E-12)
      )

type SolveFun wrapper m
  = wrapper -- ^ S
 -> V (PrimState m) ℂ -- ^ b
 -> V (PrimState m) ℂ -- ^ x
 -> m SolverStats

type MakeSFun wrapper m = Int -> M 'Row (PrimState m) ℂ -> m wrapper

-- | Sparse representation of the @S@ matrix.
data SMatrix s =
  SMatrix !(M 'Row s ℂ) -- ^ Derivatives (O - 〈O〉)
          !(V s ℂ) -- ^ Workspace of size #steps.
          !(Maybe ℂ) -- ^ Regulariser λ

makeS :: PrimMonad m => Maybe (Int -> ℂ) -> MakeSFun (SMatrix (PrimState m)) m
makeS regulariser i gradients = do
  workspace <- newDenseVector (gradients ^. dim ^. _1)
  return $! SMatrix gradients workspace ((\f -> f i) <$> regulariser)

makeS' :: PrimMonad m => Maybe (Int -> ℂ) -> MakeSFun (M 'Row (PrimState m) ℂ) m
makeS' regulariser i gradients = do
  let (steps, params) = gradients ^. dim
  s <- newDenseMatrix params params
  herk MatUpper ConjTranspose (1 / int2Float steps) gradients 0 s
  herk MatLower ConjTranspose (1 / int2Float steps) gradients 0 s
  case regulariser of
    Just f  -> do
      one <- MDenseVector @'Direct params 0 <$> V.unsafeThaw (V.singleton 1)
      axpy (f i) one (MDenseVector @'Direct params (s ^. stride + 1) (s ^. buffer))
    Nothing -> return ()
  return $! s

opS :: PrimMonad m => SMatrix (PrimState m) -> Operator m ℂ
opS (SMatrix o temp λ) x out = do
  let scale = (1 / int2Float (o ^. dim ^. _1)) :+ 0
  gemv NoTranspose 1 o x 0 temp
  case λ of
    (Just λ') -> copy x out >> gemv ConjTranspose scale o temp λ' out
    Nothing   -> gemv ConjTranspose scale o temp 0 out

{-
solveS :: PrimMonad m => CGConfig ℝ -> SolveFun (MDenseMatrix 'Row (PrimState m) ℂ) m
solveS (CGConfig maxIter tol) = doSolve
  where doSolve s b x = do
          let operator = \input output -> gemv NoTranspose 1 s input 0 output
          zero <- MDenseVector @'Direct (b ^. dim) 0 <$> V.unsafeThaw (V.singleton 0)
          copy zero x
          !answer <- CG.cg maxIter tol operator b x
          trace (show answer) (return ())
-}

solveS :: CGConfig ℝ -> SolveFun (SMatrix RealWorld) IO
solveS (CGConfig maxIter tol) = doSolve
 where
  doSolve s b x = do
    t1  <- getCPUTime
    one <- MDenseVector @ 'Direct (b ^. dim) 0 <$> V.unsafeThaw (V.singleton 0)
    copy one x
    (iters, err) <- CG.cg maxIter tol (opS s) b x
    t2           <- getCPUTime
    return $! SolverStats iters err (fromIntegral (t2 - t1) * 1.0E-12)

solveS' :: SolveFun (MDenseMatrix 'Row RealWorld ℂ) IO
solveS' s b x = assert (x ^. stride == 1) $ do
  t1 <- getCPUTime
  copy b x
  cgelsd s (MDenseMatrix (b ^. dim) 1 1 (x ^. buffer)) (-1.0)
  t2 <- getCPUTime
  return $! SolverStats (-1) (-1.0) (fromIntegral (t2 - t1) * 1.0E-12)

type Stepper v m a = v (PrimState m) a -> m ()

numberSteps :: (Int, Int, Int) -> Int
numberSteps (low, high, step) = (high - low - 1) `div` step + 1

newWorkspace :: PrimMonad m => Int -> Int -> m (DenseWorkspace (PrimState m) ℂ)
newWorkspace nParams nSteps =
  DenseWorkspace
    <$> newDenseVector nParams
    <*> newDenseMatrix nSteps nParams
    <*> newDenseVector nParams

sr
  :: SRConfig ℂ
  -> Hamiltonian
  -> MRbm (PrimState IO)
  -> (IterInfo -> IO ())
  -> IO ()
sr config hamiltonian ψ process = newWorkspace nParams nSteps
  >>= \w -> loop w 0
 where
  nParams  = size ψ
  nSteps   = config ^. mc ^. runs * numberSteps (config ^. mc ^. steps)
  doAxpy   = axpyRbm
  doMakeS  = makeS (config ^. regulariser)
  doSolveS = solveS (config ^. cg)
  doSample = sampleRbm (config ^. mc) hamiltonian
  loop !w@(DenseWorkspace f grad δ) !i
    | i >= config ^. maxIter = return ()
    | otherwise = do
      mcStats <- doSample ψ f grad
      fNorm   <- nrm2 f
      cgStats <- doMakeS i grad >>= \s -> doSolveS s f δ
      δNorm   <- nrm2 δ
      unsafeFreezeRbm ψ
        >>= \ψ' ->
              let stats = IterInfo i ψ' mcStats cgStats fNorm δNorm
              in  process stats
      doAxpy (-(config ^. rate) i) δ ψ
      loop w (i + 1)

max_stars_repo_name: twesterhout/tcm-swarm

max_issues_repo_name: twesterhout/tcm-swarm

max_forks_repo_name: twesterhout/tcm-swarm


text:
module Models where

import Layer
import Numeric.LinearAlgebra
import Data.List
import Loss

data Model a = Model [a]

addLayer :: Model a -> a -> Model a
addLayer (Model []) layer = Model [layer]
addLayer (Model l) layer = Model (l++[layer])

genParamsList :: Model VLayer -> [Matrix R]
genParamsList (Model []) = []
genParamsList m =
  let genModelParams' res (Model []) = reverse res
      genModelParams' res (Model (mx:mxs)) =
        genModelParams' ((weight mx):res) (Model mxs) in
  genModelParams' [] m

-- update the model's params from the params list
updateModelParams :: Model VLayer -> [Matrix R] -> Model VLayer
updateModelParams (Model []) _ = Model []
updateModelParams m p  =
  let updateModelParams' (Model res) (Model []) _ = Model (reverse res)
      updateModelParams' (Model res) (Model (mx:mxs)) (px:pxs) =
        let newParams = mx {weight=px} in
        updateModelParams' (Model (newParams:res)) (Model mxs) pxs
  in updateModelParams' (Model []) m p

genGParamsList :: Model VLayer -> [Matrix R]
genGParamsList (Model []) = []
genGParamsList l =
  let genModelGParams' res (Model []) = reverse res
      genModelGParams' res (Model (mx:mxs)) =
        let gp = ( (rowNum mx)  >< (colNum mx)) [0..]::(Matrix R) in
        genModelGParams' (gp:res) (Model mxs) in
  genModelGParams' [] l

--updateModelGParams :: Model VLayer -> [Matrix R] -> Model VLayer
-- update the gparams list from the model
updateGParamsList :: Model VLayer -> [Matrix R] -> [Matrix R]
updateGParamsList (Model []) _ = []
updateGParamsList l gpl =
  let updateModelGParams' res (Model []) _ = reverse res
      updateModelGParams res (Model (mx:mxs)) (gpx:gpxs) =
        let deltaUp = delta mx in
        updateModelGParams' ((deltaUp+gpx):res) (Model mxs) gpxs in
  updateModelGParams' [] l gpl

forwardInterface :: Matrix R -> VLayer -> Matrix R
forwardInterface input layer =
  case (layerType layer) of
    LinearLayer -> forwardLinearLayer input layer
    SigmoidLayer -> forwardSigmoidLayer input layer
    --SoftmaxLayer -> forwardSoftmaxLayer input layer
    ReluLayer -> forwardReluLayer input layer

backwardInterface :: VLayer -> Matrix R -> Matrix R
backwardInterface layer output =
  case (layerType layer) of
    LinearLayer -> backwardLinearLayer output layer
    SigmoidLayer -> backwardSigmoidLayer output layer
    ReluLayer -> backwardReluLayer output layer

-- forward process of the neural network model
forward :: Matrix R -> Model VLayer -> Matrix R
forward input (Model []) = input
forward input (Model xs) = foldl forwardInterface input xs

-- the backpropagate algorithm of Neural Network
-- update the delta of each layer from end to start
backward :: String -> Matrix R -> Matrix R -> Model VLayer -> Model VLayer
backward _ _ _ (Model []) = Model []
backward lossType yMat tMat l =
    let loss = getLoss lossType yMat tMat in
    let backward' (Model res) (Model []) _ = Model (reverse res)
        backward' (Model res) (Model (mx:mxs)) dw =
          let dwUp = backwardInterface mx dw
              layerUp = mx {delta=dwUp} in
          backward' (Model (layerUp:res)) (Model mxs) dwUp in
    backward' (Model []) l loss

max_stars_repo_name: neutronest/vortex

max_issues_repo_name: neutronest/vortex

max_forks_repo_name: neutronest/vortex


text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE KindSignatures        #-}
{-# LANGUAGE LambdaCase            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE RoleAnnotations       #-}
{-# LANGUAGE DuplicateRecordFields #-}

{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

--TODO figure out what export
module ObjRSimplex (
  VSimplex
  , ObjRF(..)
  , ObjRVtxOf
  , ObjRVtx(..)
  , FObjRVtx(..)
  , ObjRSimplex(..)
  , SortedObjRSimplex(..)
  , CmpObjRVtx
  , makeFObjRVtx
  , cmpObjRVtx
  , RVtxObjCmp
  , getRVtxObjCmp
  ) where

import Data.Ord (comparing)
--import Data.Function (on)
--import Data.Bool (bool)
--import Data.Maybe (fromMaybe)
--import qualified Data.Foldable as FLD
import Data.Coerce

import qualified Data.Finite as F

import qualified Data.Vector.Sized as VS

import Numeric.LinearAlgebra.Static

import GHC.TypeNats

import qualified GDP

import Data.Vector.GDP.Sorted

--TODO need another file that has SortedSimplex stuff
--must create sorted simplex with compartor that is derived from objective
--TODO problem, compartor is completely determined by ObjF in this specific case

--note this has all the axioms of Objective Simplex / Vertex, others should derive from this
--rename GDP, rename to RnObjectiveSimplex

type VSimplex v n = VS.Vector (n + 1) (v n)

type ObjRF n = R n -> ℝ

newtype ObjRVtxOf fObjName name = ObjVtxOf GDP.Defn
type role ObjRVtxOf nominal nominal

data ObjRVtx n = ObjRVtx { getX :: !(R n)
                         , getFx :: ℝ
                         }
                         deriving Show

newtype FObjRVtx fObjName n =
  FObjRVtx { getObjRVtx :: ObjRVtx n GDP.? ObjRVtxOf fObjName }

newtype ObjRSimplex fObjName n =
  ObjRSimplex { vtxs :: VSimplex (FObjRVtx fObjName) n }

newtype SortedObjRSimplex fObjName n =
  SortedObjRSimplex { vtxs :: VSimplex (FObjRVtx fObjName) n GDP.? SortedBy (RVtxObjCmp fObjName) }

type CmpObjRVtx fObjName n = Comparison (FObjRVtx fObjName n)

makeFObjRVtx :: (ObjRF n GDP.~~ fObjName) -> R n -> FObjRVtx fObjName n
makeFObjRVtx (GDP.The f) x = FObjRVtx $ GDP.assert ObjRVtx { getX = x, getFx = f x }

--TODO rewrite using GDP stripping proof function
--TODO note objective will be in environment, make new copy of nelderMead
cmpObjRVtx :: CmpObjRVtx fObjName n
cmpObjRVtx (FObjRVtx (GDP.The vtx1)) (FObjRVtx (GDP.The vtx2)) =
  comparing getFx vtx1 vtx2

newtype RVtxObjCmp fObjName = RVtxObjCmp GDP.Defn
type role RVtxObjCmp nominal

getRVtxObjCmp :: (ObjRF n GDP.~~ fObjName)
  -> (CmpObjRVtx fObjName n GDP.~~ RVtxObjCmp fObjName)
getRVtxObjCmp f = GDP.defn cmpObjRVtx

newtype VtxsOfObjRSimplex vName = VtxsOfObjRSimplex GDP.Defn
type role VtxsOfObjRSimplex nominal

getObjRSimplexNamedVtxs :: forall fObjName smplxName n.
  ObjRSimplex fObjName n GDP.~~ smplxName
  -> (VSimplex (FObjRVtx fObjName) n GDP.~~ VtxsOfObjRSimplex smplxName)
getObjRSimplexNamedVtxs (GDP.The s) = GDP.defn $ vtxs (s :: ObjRSimplex fObjName n)

max_stars_repo_name: brnzhg/howod

max_issues_repo_name: brnzhg/howod

max_forks_repo_name: brnzhg/howod


text:
{-# LANGUAGE OverloadedStrings, NoImplicitPrelude, BangPatterns #-}
module Main where
import ClassyPrelude
import qualified Numeric.LinearAlgebra as HMatrix
import NLP.Albemarle.Dict (counts, ids)
import qualified NLP.Albemarle.Dict as Dict
import NLP.Albemarle.LSA (termvectors, topicweights)
import qualified NLP.Albemarle.LSA as LSA
import qualified NLP.Albemarle.Tokens as Tokens
import qualified System.IO.Streams as Streams
import Data.Text.Encoding.Error (lenientDecode)
import Lens.Micro

main :: IO ()
main = do
  (final_dict, model) <- Streams.withFileAsInput "kjv.verses.txt.gz" $ \file ->
    Streams.gunzip file
    >>= Streams.lines
    >>= Streams.decodeUtf8With lenientDecode
    >>= Streams.map Tokens.wordTokenize
    >>= Streams.chunkList 5000 -- Yep, I'm making this stuff up.
    >>= Streams.map (\chunk -> let
      dict = Dict.dictifyAllWords chunk
      sparsem = Dict.asSparseMatrix dict chunk
      lsa = LSA.lsa 100 sparsem
      in (dict, lsa))
    >>= Streams.fold (\(!d1, !lsa1) (d2, lsa2) -> let
      d3 = Dict.filterDict 2 0.5 1000 $ d1 <> d2
      lsa3 = LSA.rebase d1 d3 lsa1 <> LSA.rebase d2 d3 lsa2
      in (d3, lsa3)) (mempty, mempty)

  -- It should use all 100 words allowed plus the unknown
  print $ length (final_dict^.counts)
  -- It should have a full size LSA as well
  print $ HMatrix.size (model^.termvectors)
  -- Plus topic weights
  print $ HMatrix.size (model^.topicweights)

max_stars_repo_name: SeanTater/albemarle

max_issues_repo_name: SeanTater/albemarle

max_forks_repo_name: SeanTater/albemarle


text:
{-
   Copyright 2016, Dominic Orchard, Andrew Rice, Mistral Contrastin, Matthew Danish

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
-}

{-
  Units of measure extension to Fortran: backend
-}

{-# LANGUAGE TupleSections #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Camfort.Specification.Units.InferenceBackend
  ( chooseImplicitNames
  , criticalVariables
  , inconsistentConstraints
  , inferVariables
  -- mainly for debugging and testing:
  , shiftTerms
  , flattenConstraints
  , flattenUnits
  , constraintsToMatrix
  , constraintsToMatrices
  , rref
  , genUnitAssignments
  , genUnitAssignments'
  , provenance
  , splitNormHNF
  ) where

import           Camfort.Specification.Units.Environment
import qualified Camfort.Specification.Units.InferenceBackendFlint as Flint
import           Control.Arrow (first, second, (***))
import           Control.Monad
import           Control.Monad.ST
import           Control.Parallel.Strategies
import qualified Data.Array as A
import           Data.Generics.Uniplate.Operations
  (transformBi, universeBi)
import           Data.Graph.Inductive hiding ((><))
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import           Data.List
  ((\\), findIndex, inits, nub, partition, sort, sortBy, group, tails, foldl')
import qualified Data.Map.Strict as M
import           Data.Maybe (fromMaybe, mapMaybe)
import           Data.Ord
import           Data.Tuple (swap)
import           Numeric.LinearAlgebra
  ( atIndex, (<>)
  , rank, (?), (¿)
  , rows, cols
  , subMatrix, diag
  , fromBlocks, ident
  )
import qualified Numeric.LinearAlgebra as H
import           Numeric.LinearAlgebra.Devel
  ( newMatrix, readMatrix
  , writeMatrix, runSTMatrix
  , freezeMatrix, STMatrix
  )
import           Prelude hiding ((<>))

-- | Returns list of formerly-undetermined variables and their units.
inferVariables :: Constraints -> [(VV, UnitInfo)]
inferVariables cons = unitVarAssignments
  where
    unitAssignments = genUnitAssignments cons
    -- Find the rows corresponding to the distilled "unit :: var"
    -- information for ordinary (non-polymorphic) variables.
    unitVarAssignments            =
      [ (var, units) | ([UnitPow (UnitVar var)                 k], units) <- unitAssignments, k `approxEq` 1 ] ++
      [ (var, units) | ([UnitPow (UnitParamVarAbs (_, var)) k], units)    <- unitAssignments, k `approxEq` 1 ]

-- Detect inconsistency if concrete units are assigned an implicit
-- abstract unit variable with coefficients not equal, or there are
-- monomorphic literals being given parametric polymorphic units.
detectInconsistency :: [([UnitInfo], UnitInfo)] -> Constraints
detectInconsistency unitAssignments = inconsist
  where
    ua' = map (shiftTerms . fmap flattenUnits) unitAssignments
    badImplicits = [ fmap foldUnits a | a@([UnitPow (UnitParamImpAbs _) k1], rhs) <- ua'
                                      , UnitPow _ k2 <- rhs
                                      , k1 /= k2 ]
    inconsist = unitAssignmentsToConstraints badImplicits ++ mustBeUnitless unitAssignments

-- Must be unitless: any assignments of parametric abstract units to
-- monomorphic literals.
mustBeUnitless :: [([UnitInfo], UnitInfo)] -> Constraints
mustBeUnitless unitAssignments = mbu
  where
    mbu = [ ConEq UnitlessLit (UnitPow (UnitLiteral l) k)
          | (UnitPow (UnitLiteral l) k:_, rhs) <- ua''
          , any isParametric (universeBi rhs :: [UnitInfo]) ]
    -- ua' = map (shiftTerms . fmap flattenUnits) unitAssignments
    ua'' = map (shiftTermsBy isLiteral . fmap flattenUnits) unitAssignments

    isLiteral UnitLiteral{} = True
    isLiteral (UnitPow UnitLiteral{} _) = True
    isLiteral _ = False

    isParametric UnitParamVarAbs{} = True
    isParametric UnitParamPosAbs{} = True
    isParametric UnitParamEAPAbs{} = True
    isParametric UnitParamLitAbs{} = True
    isParametric UnitParamImpAbs{} = True
    isParametric (UnitPow u _)     = isParametric u
    isParametric _                 = False

-- convert the assignment format back into constraints
unitAssignmentsToConstraints :: [([UnitInfo], UnitInfo)] -> Constraints
unitAssignmentsToConstraints = map (uncurry ConEq . first foldUnits)

-- | Raw units-assignment pairs.
genUnitAssignments :: Constraints -> [([UnitInfo], UnitInfo)]
genUnitAssignments cons
  -- if the results include any mappings that must be forced to be unitless...
  | mbu <- mustBeUnitless ua, not (null mbu) = genUnitAssignments (mbu ++ unitAssignmentsToConstraints ua)
  | null (detectInconsistency ua)            = ua
  | otherwise                                = []
  where
    ua = genUnitAssignments' colSort cons

-- | Break up the problem of solving normHNF on each group of related
-- columns, then bring it all back together.
splitNormHNF :: H.Matrix Double -> (H.Matrix Double, [Int])
splitNormHNF unsolvedM = (combinedMat, allNewColIndices)
  where
    combinedMat      = joinMat (map (first fst) solvedMs)
    allNewColIndices = concatMap (snd . fst) solvedMs

    inParallel = (`using` parTuple2 (parList rseq) rseq)
    (solvedMs, _) = inParallel . foldl' eachResult ([], cols unsolvedM) $ map (first Flint.normHNF) (splitMat unsolvedM)

    -- for each result re-number the generated columns & add mappings for each.
    eachResult (ms, startI) ((m, newColIndices), origCols) = (((m, newColIndices'), origCols'):ms, endI)
      where
        -- produce (length newColIndices) number of mappings
        endI           = startI + length newColIndices
        -- re-number the newColIndices according to the lookup list
        newColIndices' = map (origCols !!) newColIndices
        -- add columns in the (combined) matrix for the newly
        -- generated columns from running normHNF on m.
        origCols'      = origCols ++ [startI .. endI-1]

genUnitAssignments' :: SortFn -> Constraints -> [([UnitInfo], UnitInfo)]
genUnitAssignments' _ [] = []
genUnitAssignments' sortfn cons
  | null colList                                      = []
  | null inconsists                                   = unitAssignments
  | otherwise                                         = []
  where
    (lhsM, rhsM, inconsists, lhsColA, rhsColA) = constraintsToMatrices' sortfn cons
    unsolvedM | rows rhsM == 0 || cols rhsM == 0 = lhsM
              | rows lhsM == 0 || cols lhsM == 0 = rhsM
              | otherwise                        = fromBlocks [[lhsM, rhsM]]
    (solvedM, newColIndices)      = splitNormHNF unsolvedM
    -- solvedM can have additional columns and rows from normHNF;
    -- cosolvedM corresponds to the original lhsM.
    -- cosolvedM                     = subMatrix (0, 0) (rows solvedM, cols lhsM) solvedM
    -- cosolvedMrhs                  = subMatrix (0, cols lhsM) (rows solvedM, cols solvedM - cols lhsM) solvedM

    -- generate a colList with both the original columns and new ones generated
    -- if a new column generated was derived from the right-hand side then negate it
    numLhsCols                    = 1 + snd (A.bounds lhsColA)
    colList                       = map (1,) (A.elems lhsColA ++ A.elems rhsColA) ++ map genC newColIndices
    genC n | n >= numLhsCols      = (-k, UnitParamImpAbs (show u))
           | otherwise            = (k, UnitParamImpAbs (show u))
      where (k, u) = colList !! n
    -- Convert the rows of the solved matrix into flattened unit
    -- expressions in the form of "unit ** k".
    unitPow (k, u) x              = UnitPow u (k * x)
    unitPows                      = map (concatMap flattenUnits . zipWith unitPow colList) (H.toLists solvedM)

    -- Variables to the left, unit names to the right side of the equation.
    unitAssignments               = map (fmap (foldUnits . map negatePosAbs) . checkSanity . partition (not . isUnitRHS')) unitPows
    isUnitRHS' (UnitPow (UnitName _) _)        = True
    isUnitRHS' (UnitPow (UnitParamEAPAbs _) _) = True
    -- Because this version of isUnitRHS different from
    -- constraintsToMatrix interpretation, we need to ensure that any
    -- moved ParamPosAbs units are negated, because they are
    -- effectively being shifted across the equal-sign:
    isUnitRHS' (UnitPow (UnitParamImpAbs _) _) = True
    isUnitRHS' (UnitPow (UnitParamPosAbs (_, 0)) _) = False
    isUnitRHS' (UnitPow (UnitParamPosAbs _) _) = True
    isUnitRHS' _                               = False

checkSanity :: ([UnitInfo], [UnitInfo]) -> ([UnitInfo], [UnitInfo])
checkSanity (u1@[UnitPow (UnitVar _) _], u2)
  | or $ [ True | UnitParamPosAbs (_, _) <- universeBi u2 ]
      ++ [ True | UnitParamImpAbs _      <- universeBi u2 ] = (u1++u2,[])
checkSanity (u1@[UnitPow (UnitParamVarAbs (f, _)) _], u2)
  | or [ True | UnitParamPosAbs (f', _) <- universeBi u2, f' /= f ] = (u1++u2,[])
checkSanity c = c

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

-- FIXME: you know better...
approxEq :: Double -> Double -> Bool
approxEq a b = abs (b - a) < epsilon
notApproxEq :: Double -> Double -> Bool
notApproxEq a b = not (approxEq a b)
epsilon :: Double
epsilon = 0.001 -- arbitrary

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

type RowNum = Int               -- ^ 'row number' of matrix
type ColNum = Int               -- ^ 'column number' of matrix
-- | Represents a subproblem of AX=B where the row numbers and column
-- numbers help you re-map back to the original problem.
type Subproblem = ([RowNum], (H.Matrix Double, H.Matrix Double), [ColNum])

-- | Divide up the AX=B problem into smaller problems based on the
-- 'related columns' and their corresponding rows in the
-- right-hand-side of the equation. Where lhsM = A and rhsM = B.  The
-- resulting list of subproblems contains the new, smaller As and Bs
-- as well as a list of original row numbers and column numbers to
-- aide re-mapping back to the original lhsM and rhsM.
splitMatWithRHS :: H.Matrix Double -> H.Matrix Double -> [Subproblem]
splitMatWithRHS lhsM rhsM | cols lhsM > 0 = map (eachComponent . sort) $ scc (relatedColumnsGraph lhsM)
                          | otherwise     = []
  where
    -- Gets called on every strongly-connected component / related set of columns.
    eachComponent cs = (rs, mats, cs)
      where
        -- Selected columns
        lhsSelCols :: H.Matrix Double
        lhsSelCols = lhsM ¿ cs

        csLen = cols lhsSelCols

        -- Find the row numbers of the 'all zero' rows in lhsM.
        lhsAllZeroRows :: [RowNum]
        lhsAllZeroRows = map fst . filter (all (approxEq 0) . snd) . zip [0..] $ H.toLists lhsM

        -- Find the row numbers that correspond to the non-zero co-efficients in the selected columns.
        lhsNonZeroColRows :: [(RowNum, [Double])]
        lhsNonZeroColRows = filter (any (notApproxEq 0) . snd) . zip [0..] . H.toLists $ lhsSelCols

        -- List of all the row numbers and row values combined from the two above variables.
        lhsNumberedRows :: [(RowNum, [Double])]
        lhsNumberedRows = sortBy (comparing fst) $ lhsNonZeroColRows ++ zip lhsAllZeroRows (repeat (replicate csLen 0))

        -- For each of the above LHS rows find a corresponding RHS row.
        rhsSelRows :: [[Double]]
        rhsSelRows | rows rhsM > 0 = H.toLists (rhsM ? map fst lhsNumberedRows)
                   | otherwise     = []

        reassoc (a, b) c = (a, (b, c))

        notAllZero (_, (lhs, rhs)) = any (notApproxEq 0) (lhs ++ rhs)

        -- Zip the selected LHS, RHS rows together, filter out any that are all zeroes.
        numberedRows :: ([RowNum], [([Double], [Double])])
        numberedRows = unzip . filter notAllZero $ zipWith reassoc lhsNumberedRows rhsSelRows

        rs :: [RowNum]          -- list of row numbers in the subproblem
        mats :: (H.Matrix Double, H.Matrix Double) -- LHS/RHS subproblem matrices
        (rs, mats) = second ((H.fromLists *** H.fromLists) . unzip) numberedRows

-- | Split the lhsM/rhsM problem into subproblems and then look for
-- inconsistent rows in each subproblem, concatenating all of the
-- inconsistent row numbers found (in terms of the rows of the
-- original lhsM).
splitFindInconsistentRows :: H.Matrix Double -> H.Matrix Double -> [RowNum]
splitFindInconsistentRows lhsMat rhsMat = concatMap eachComponent $ splitMatWithRHS lhsMat rhsMat
  where
    eachComponent (rs, (lhsM, rhsM), _) = map (rs !!) $ findInconsistentRows lhsM augM
      where
        -- Augmented matrix is defined as the combined LHS/RHS matrices.
        augM
          | rows rhsM == 0 || cols rhsM == 0 = lhsM
          | rows lhsM == 0 || cols lhsM == 0 = rhsM
          | otherwise = fromBlocks [[lhsM, rhsM]]

-- | Break out the 'unrelated' columns in a single matrix into
-- separate matrices, along with a list of their original column
-- positions.
splitMat :: H.Matrix Double -> [(H.Matrix Double, [ColNum])]
splitMat m = map (eachComponent . sort) $ scc (relatedColumnsGraph m)
  where
    eachComponent cs = (H.fromLists . filter (any (/= 0)) . H.toLists $ m ¿ cs, cs)

-- | Bring together the split matrices and put the columns back in
-- their original order. Rows may not be in the same order as the
-- original, but the constraints should be equivalent.
joinMat :: [(H.Matrix Double, [Int])] -> H.Matrix Double
joinMat ms = sortedM
  where
    disorderedM = H.diagBlock (map fst ms)
    colsWithIdx = zip (concatMap snd ms) . H.toColumns $ disorderedM
    sortedM     = H.fromColumns . map snd . sortBy (comparing fst) $ colsWithIdx

-- | Turn a matrix into a graph where each node represents a column
-- and each edge represents two columns that have non-zero
-- co-efficients in some row. Basically, 'related columns'. Also
-- includes self-refs for each node..
relatedColumnsGraph :: H.Matrix Double -> Gr () ()
relatedColumnsGraph m = mkGraph (map (,()) ns) (map (\ (a,b) -> (a,b,())) es)
  where
    nonZeroCols = [ [ j | j <- [0..cols m - 1], not (m `atIndex` (i, j) `approxEq` 0) ] | i <- [0..rows m - 1] ]
    ns          = nub $ concat nonZeroCols
    es          = [ (i, j) | cs <- nonZeroCols, [i, j] <- sequence [cs, cs] ]

-- Convert a set of constraints into a matrix of co-efficients, and a
-- reverse mapping of column numbers to units.
constraintsToMatrix :: Constraints -> (H.Matrix Double, [Int], A.Array Int UnitInfo)
constraintsToMatrix cons
  | all null lhs = (H.ident 0, [], A.listArray (0, -1) [])
  | otherwise = (augM, inconsists, A.listArray (0, length colElems - 1) colElems)
  where
    -- convert each constraint into the form (lhs, rhs)
    consPairs       = filter (uncurry (/=)) $ flattenConstraints cons
    -- ensure terms are on the correct side of the equal sign
    shiftedCons     = map shiftTerms consPairs
    lhs             = map fst shiftedCons
    rhs             = map snd shiftedCons
    (lhsM, lhsCols) = flattenedToMatrix colSort lhs
    (rhsM, rhsCols) = flattenedToMatrix colSort rhs
    colElems        = A.elems lhsCols ++ A.elems rhsCols
    augM            = if rows rhsM == 0 || cols rhsM == 0 then lhsM else if rows lhsM == 0 || cols lhsM == 0 then rhsM else fromBlocks [[lhsM, rhsM]]
    inconsists      = splitFindInconsistentRows lhsM rhsM

constraintsToMatrices :: Constraints -> (H.Matrix Double, H.Matrix Double, [Int], A.Array Int UnitInfo, A.Array Int UnitInfo)
constraintsToMatrices cons = constraintsToMatrices' colSort cons

constraintsToMatrices' :: SortFn -> Constraints -> (H.Matrix Double, H.Matrix Double, [Int], A.Array Int UnitInfo, A.Array Int UnitInfo)
constraintsToMatrices' sortfn cons
  | all null lhs = (H.ident 0, H.ident 0, [], A.listArray (0, -1) [], A.listArray (0, -1) [])
  | otherwise = (lhsM, rhsM, inconsists, lhsCols, rhsCols)
  where
    -- convert each constraint into the form (lhs, rhs)
    consPairs       = filter (uncurry (/=)) $ flattenConstraints cons
    -- ensure terms are on the correct side of the equal sign
    shiftedCons     = map shiftTerms consPairs
    lhs             = map fst shiftedCons
    rhs             = map snd shiftedCons
    (lhsM, lhsCols) = flattenedToMatrix sortfn lhs
    (rhsM, rhsCols) = flattenedToMatrix sortfn rhs
    inconsists      = splitFindInconsistentRows lhsM rhsM

-- [[UnitInfo]] is a list of flattened constraints
flattenedToMatrix :: SortFn -> [[UnitInfo]] -> (H.Matrix Double, A.Array Int UnitInfo)
flattenedToMatrix sortfn cons = (m, A.array (0, numCols - 1) (map swap uniqUnits))
  where
    m = runSTMatrix $ do
          newM <- newMatrix 0 numRows numCols
          -- loop through all constraints
          forM_ (zip cons [0..]) $ \ (unitPows, row) -> do
            -- write co-efficients for the lhs of the constraint
            forM_ unitPows $ \ (UnitPow u k) -> do
              case M.lookup u colMap of
                Just col -> readMatrix newM row col >>= (writeMatrix newM row col . (+k))
                _        -> return ()
          return newM
    -- identify and enumerate every unit uniquely
    uniqUnits = flip zip [0..] . map head . group . sortBy sortfn $ [ u | UnitPow u _ <- concat cons ]
    -- map units to their unique column number
    colMap    = M.fromList uniqUnits
    numRows   = length cons
    numCols   = M.size colMap

negateCons :: [UnitInfo] -> [UnitInfo]
negateCons = map (\ (UnitPow u k) -> UnitPow u (-k))

negatePosAbs :: UnitInfo -> UnitInfo
negatePosAbs (UnitPow (UnitParamPosAbs x) k) = UnitPow (UnitParamPosAbs x) (-k)
negatePosAbs (UnitPow (UnitParamImpAbs v) k) = UnitPow (UnitParamImpAbs v) (-k)
negatePosAbs u                               = u

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

-- Units that should appear on the right-hand-side of the matrix during solving
isUnitRHS :: UnitInfo -> Bool
isUnitRHS (UnitPow (UnitName _) _)        = True
isUnitRHS (UnitPow (UnitParamEAPAbs _) _) = True
isUnitRHS _                               = False

-- | Shift UnitNames/EAPAbs poly units to the RHS, and all else to the LHS.
shiftTerms :: ([UnitInfo], [UnitInfo]) -> ([UnitInfo], [UnitInfo])
shiftTerms (lhs, rhs) = (lhsOk ++ negateCons rhsShift, rhsOk ++ negateCons lhsShift)
  where
    (lhsOk, lhsShift) = partition (not . isUnitRHS) lhs
    (rhsOk, rhsShift) = partition isUnitRHS rhs

-- | Shift terms based on function f (<- True, False ->).
shiftTermsBy :: (UnitInfo -> Bool) -> ([UnitInfo], [UnitInfo]) -> ([UnitInfo], [UnitInfo])
shiftTermsBy f (lhs, rhs) = (lhsOk ++ negateCons rhsShift, rhsOk ++ negateCons lhsShift)
  where
    (lhsOk, lhsShift) = partition f lhs
    (rhsOk, rhsShift) = partition (not . f) rhs

-- | Translate all constraints into a LHS, RHS side of units.
flattenConstraints :: Constraints -> [([UnitInfo], [UnitInfo])]
flattenConstraints = map (\ (ConEq u1 u2) -> (flattenUnits u1, flattenUnits u2))

--------------------------------------------------
-- Matrix solving functions based on HMatrix

-- | Returns given matrix transformed into Reduced Row Echelon Form
rref :: H.Matrix Double -> H.Matrix Double
rref a = snd $ rrefMatrices' a 0 0 []
  where
    -- (a', den, r) = Flint.rref a

-- Provenance of matrices.
data RRefOp
  = ElemRowSwap Int Int         -- ^ swapped row with row
  | ElemRowMult Int Double      -- ^ scaled row by constant
  | ElemRowAdds [(Int, Int)]    -- ^ set of added row onto row ops
  deriving (Show, Eq, Ord)

-- worker function
-- invariant: the matrix a is in rref except within the submatrix (j-k,j) to (n,n)
rrefMatrices' :: H.Matrix Double -> Int -> Int -> [(H.Matrix Double, RRefOp)] ->
                 ([(H.Matrix Double, RRefOp)], H.Matrix Double)
rrefMatrices' a j k mats
  -- Base cases:
  | j - k == n            = (mats, a)
  | j     == m            = (mats, a)

  -- When we haven't yet found the first non-zero number in the row, but we really need one:
  | a `atIndex` (j - k, j) == 0 = case findIndex (/= 0) below of
    -- this column is all 0s below current row, must move onto the next column
    Nothing -> rrefMatrices' a (j + 1) (k + 1) mats
    -- we've found a row that has a non-zero element that can be swapped into this row
    Just i' -> rrefMatrices' (swapMat <> a) j k ((swapMat, ElemRowSwap i (j - k)):mats)
      where i       = j - k + i'
            swapMat = elemRowSwap n i (j - k)

  -- We have found a non-zero cell at (j - k, j), so transform it into
  -- a 1 if needed using elemRowMult, and then clear out any lingering
  -- non-zero values that might appear in the same column, using
  -- elemRowAdd:
  | otherwise             = rrefMatrices' a2 (j + 1) k mats2
  where
    n     = rows a
    m     = cols a
    below = getColumnBelow a (j - k, j)
    scale = recip (a `atIndex` (j - k, j))
    erm   = elemRowMult n (j - k) scale

    -- scale the row if the cell is not already equal to 1
    (a1, mats1) | a `atIndex` (j - k, j) /= 1 = (erm <> a, (erm, ElemRowMult (j - k) scale):mats)
                | otherwise                   = (a, mats)

    -- Locate any non-zero values in the same column as (j - k, j) and
    -- cancel them out. Optimisation: instead of constructing a
    -- separate elemRowAdd matrix for each cancellation that are then
    -- multiplied together, simply build a single matrix that cancels
    -- all of them out at the same time, using the ST Monad.
    findAdds _ curM ms
      | isWritten = (newMat <> curM, (newMat, ElemRowAdds matOps):ms)
      | otherwise = (curM, ms)
      where
        (isWritten, matOps, newMat) = runST $ do
          newM <- newMatrix 0 n n :: ST s (STMatrix s Double)
          sequence_ [ writeMatrix newM i' i' 1 | i' <- [0 .. (n - 1)] ]
          let f w o i | i >= n                  = return (w, o)
                      | i == j - k              = f w o (i + 1)
                      | a `atIndex` (i, j) == 0 = f w o (i + 1)
                      | otherwise               = writeMatrix newM i (j - k) (- (a `atIndex` (i, j)))
                                                  >> f True ((i, j - k):o) (i + 1)
          (isW, ops) <- f False [] 0
          (isW, ops,) `fmap` freezeMatrix newM

    (a2, mats2) = findAdds (0::Int) a1 mats1

-- Get a list of values that occur below (i, j) in the matrix a.
getColumnBelow :: H.Matrix Double -> (Int, Int) -> [Double]
getColumnBelow a (i, j) = concat . H.toLists $ subMatrix (i, j) (n - i, 1) a
  where n = rows a

-- 'Elementary row operation' matrices
elemRowMult :: Int -> Int -> Double -> H.Matrix Double
elemRowMult n i k = diag (H.fromList (replicate i 1.0 ++ [k] ++ replicate (n - i - 1) 1.0))

elemRowSwap :: Int -> Int -> Int -> H.Matrix Double
elemRowSwap n i j
  | i == j          = ident n
  | i > j           = elemRowSwap n j i
  | otherwise       = ident n ? ([0..i-1] ++ [j] ++ [i+1..j-1] ++ [i] ++ [j+1..n-1])

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

type GraphCol = IM.IntMap IS.IntSet   -- graph from origin to dest.
type Provenance = IM.IntMap IS.IntSet -- graph from dest. to origin

opToGraphCol :: RRefOp -> GraphCol
opToGraphCol ElemRowMult{} = IM.empty
opToGraphCol (ElemRowSwap i j) = IM.fromList [ (i, IS.singleton j), (j, IS.singleton i) ]
opToGraphCol (ElemRowAdds l)   = IM.fromList $ concat [ [(i, IS.fromList [i,j]), (j, IS.singleton j)]  | (i, j) <- l ]

graphColCombine :: GraphCol -> GraphCol -> GraphCol
graphColCombine g1 g2 = IM.unionWith (curry snd) g1 $ IM.map (IS.fromList . trans . IS.toList) g2
  where
    trans = concatMap (\ i -> [i] `fromMaybe` (IS.toList <$> IM.lookup i g1))

invertGraphCol :: GraphCol -> GraphCol
invertGraphCol g = IM.fromListWith IS.union [ (i, IS.singleton j) | (j, jset) <- IM.toList g, i <- IS.toList jset ]

provenance :: H.Matrix Double -> (H.Matrix Double, Provenance)
provenance m = (m', p)
  where
    (matOps, m') = rrefMatrices' m 0 0 []
    p = invertGraphCol . foldl' graphColCombine IM.empty . map opToGraphCol $ map snd matOps

-- Worker functions:

findInconsistentRows :: H.Matrix Double -> H.Matrix Double -> [Int]
findInconsistentRows coA augA | rows augA < 2 = []
                              | otherwise     = inconsistent
  where
    inconsistent = [0..(rows augA - 1)] \\ consistent

    consistent
      -- if the space is relatively small, try it all
      | rows augA < 16 = head (filter tryRows (powerset $ reverse [0..(rows augA - 1)]))
      | otherwise = head (filter tryRows (tails ( [0..(rows augA - 1)])) ++ [[]])

    powerset = filterM (const [True, False])

    -- Rouché–Capelli theorem is that if the rank of the coefficient
    -- matrix is not equal to the rank of the augmented matrix then
    -- the system of linear equations is inconsistent.
    tryRows [] = True
    tryRows ns = (rank coA' == rank augA')
      where
        coA'  = coA ? ns
        augA' = augA ? ns

-- | Create unique names for all of the inferred implicit polymorphic
-- unit variables.
chooseImplicitNames :: [(VV, UnitInfo)] -> [(VV, UnitInfo)]
chooseImplicitNames vars = replaceImplicitNames (genImplicitNamesMap vars) vars

genImplicitNamesMap :: Data a => a -> M.Map UnitInfo UnitInfo
genImplicitNamesMap x = M.fromList [ (absU, UnitParamEAPAbs (newN, newN)) | (absU, newN) <- zip absUnits newNames ]
  where
    absUnits = nub [ u | u@(UnitParamPosAbs _)             <- universeBi x ] ++
               nub [ u | u@(UnitParamImpAbs _)             <- universeBi x ]
    eapNames = nub $ [ n | (UnitParamEAPAbs (_, n))      <- universeBi x ] ++
                     [ n | (UnitParamEAPUse ((_, n), _)) <- universeBi x ]
    newNames = filter (`notElem` eapNames) . map ('\'':) $ nameGen
    nameGen  = concatMap sequence . tail . inits $ repeat ['a'..'z']

replaceImplicitNames :: Data a => M.Map UnitInfo UnitInfo -> a -> a
replaceImplicitNames implicitMap = transformBi replace
  where
    replace u@(UnitParamPosAbs _) = fromMaybe u $ M.lookup u implicitMap
    replace u@(UnitParamImpAbs _) = fromMaybe u $ M.lookup u implicitMap
    replace u                     = u

-- | Identifies the variables that need to be annotated in order for
-- inference or checking to work.
criticalVariables :: Constraints -> [UnitInfo]
criticalVariables [] = []
criticalVariables cons = filter (not . isUnitRHS') $ map (colA A.!) criticalIndices
  where
    (unsolvedM, _, colA)          = constraintsToMatrix cons
    solvedM                       = rref unsolvedM
    uncriticalIndices             = mapMaybe (findIndex (/= 0)) $ H.toLists solvedM
    criticalIndices               = A.indices colA \\ uncriticalIndices
    isUnitRHS' (UnitName _)       = True; isUnitRHS' _ = False

-- | Returns just the list of constraints that were identified as
-- being possible candidates for inconsistency, if there is a problem.
inconsistentConstraints :: Constraints -> Maybe Constraints
inconsistentConstraints [] = Nothing
inconsistentConstraints cons
  | not (null direct) = Just direct
  | null inconsists   = Nothing
  | otherwise         = Just [ con | (con, i) <- zip cons [0..], i `elem` inconsists ]
  where
    (_, _, inconsists, _, _) = constraintsToMatrices cons
    direct = detectInconsistency $ genUnitAssignments' colSort cons

max_stars_repo_name: raehik/camfort

max_issues_repo_name: apthorpe/camfort

max_forks_repo_name: raehik/camfort


text:
{-# LANGUAGE DataKinds        #-}
{-# LANGUAGE TypeApplications #-}

module Lib.STFTSpec where

import           Data.Coerce
import           Data.Complex
import           Data.Proxy
import qualified Data.Vector     as V
import           GHC.TypeNats
import qualified Lib.STFT        as STFT
import qualified Lib.Util        as U
import           Test.Hspec
import           Test.QuickCheck

takeLast :: Int -> V.Vector a -> V.Vector a
takeLast n vec = V.reverse $ V.take n $ V.reverse vec

spec = do

  describe "STFT.hamming" $ do
    it "creates a hamming window of size N" $ do
      let result = V.fromList [8.000000000000002e-2 :+ 0.0,
              0.54 :+ 0.0,1.0 :+ 0.0,
              0.5400000000000001 :+ 0.0,
              8.000000000000002e-2 :+ 0.0
            ]
      STFT.runWindow STFT.hamming (STFT.MkWinsz (Proxy @5)) `shouldBe` result

  describe "STFT.zeroPhaseZeroPad" $ do
    it "peforms zero phase and zero padding" $ do
      tuple <- U.readWav "static/audio/singing-female.wav"
      let firstFive = V.fromList [(-0.553192138671875) :+ 0.0,
              (-0.55218505859375) :+ 0.0,
              (-0.556488037109375) :+ 0.0,
              (-0.58892822265625) :+ 0.0,
              (-0.610992431640625) :+ 0.0
            ]
          lastFive = V.fromList [(-0.443511962890625) :+ 0.0,
              (-0.479644775390625) :+ 0.0,
              (-0.49029541015625) :+ 0.0,
              (-0.493133544921875) :+ 0.0,
              (-0.52313232421875) :+ 0.0
            ]
          (_, audio) = U.makeSignal tuple
          result :: V.Vector (Complex Double)
          result = coerce $ STFT.zeroPhaseZeroPad (STFT.MkFFTsz (Proxy @2048)) audio
      V.take 5 result `shouldBe` firstFive
      takeLast 5 result `shouldBe` lastFive

  describe "STFT.dftAnal" $ do
    it "performs a dft of the given FFT size" $ do
      tuple <- U.readWav "static/audio/singing-female.wav"
      let firstFive = V.fromList [-12.049319770407699,
              -3.8232231987269523,
              -3.5161910098012816,
              -5.778847744891994,
              -4.450941161692169
            ]
          lastFive = V.fromList [9.085997084736617,
              12.384474935554943,
              11.667317158981483,
              20.354458129771224,
              11.951847223595225
            ]
          (_, audio) = U.makeSignal tuple
          (STFT.MkSpect (STFT.MkMag mag, _)) = STFT.dftAnal (STFT.MkFFTsz (Proxy @2048)) audio
      V.take 5 mag `shouldBe` firstFive
      takeLast 5 mag `shouldBe` lastFive

  describe "STFT.stft" $ do
    it "Performs the short time fourier transform" $ do
      (_, audio) <- U.readWav "static/audio/singing-female.wav"
      let firstFive = V.fromList [-95.55314985231192,
              -75.57780027962905,
              -72.63091749228488,
              -74.5060313038982,
              -80.35599835805851
            ]
          lastFive = V.fromList [-132.57867497788396,
              -128.76434452701437,
              -128.15973988408462,
              -131.54917075905306,
              -140.28246959896572
            ]
          (STFT.MkSpect (STFT.MkMag mag, _)) = head $ STFT.stft STFT.config (U.makeSig audio)
      V.take 5 mag `shouldBe` firstFive
      takeLast 5 mag `shouldBe` lastFive

  describe "STFT.divvy" $ do
    it "Splits the vector into " $ do
      let vec = V.replicate 272243 ()
          divvyVec = STFT.divvy (STFT.MkWinsz (Proxy @2048)) (STFT.MkHopsz (Proxy @256)) vec
      fmap length divvyVec `shouldBe` replicate 1056 2048

max_stars_repo_name: davlum/haskell-twm

max_issues_repo_name: davlum/haskell-twm

max_forks_repo_name: davlum/haskell-twm


text:
{-# OPTIONS_GHC -fno-warn-unticked-promoted-constructors #-}

module Noether.Algebra.Linear.Module where

import           Data.Complex

import           Noether.Lemmata.Prelude
import           Noether.Lemmata.TypeFu

import           Noether.Algebra.Actions
import           Noether.Algebra.Multiple
import           Noether.Algebra.Single

import           Noether.Algebra.Tags

data LeftModuleE
  = LeftModule_Ring_AbelianGroup_Linear_Compatible
    { leftModule_ring           :: RingE
    , leftModule_abelianGroup   :: AbelianGroupE
    , leftModule_actorLinearity :: ActorLinearE
    , leftModule_acteeLinearity :: ActeeLinearE
    , leftModule_compatibility  :: CompatibleE}
  | LeftModule_Named Symbol
                     LeftModuleE

-- | A left module (v, a) over the ring (r, p, m).
class LeftModuleK op p m r a v s

type family LeftModuleS (op :: k0) (p :: k1) (m :: k2) r (a :: k3) v :: LeftModuleE

type LeftModuleC op p m r a v = LeftModuleK op p m r a v (LeftModuleS op p m r a v)

instance ( RingK p m r zr
         , AbelianGroupK a v zag
         , ActorLinearK L m p r a v zor
         , ActeeLinearK L m r a v zee
         , CompatibleK L op m r v zlc
         ) =>
         LeftModuleK op p m r a v
           (LeftModule_Ring_AbelianGroup_Linear_Compatible zr zag zor zee zlc)

type LeftModule op p m r a v =
  ( LeftModuleC op p m r a v
  , Ring p m r
  , AbelianGroup a v
  , LinearActsOn L m p r a v
  , LeftCompatible op m r v
  )

data RightModuleE
  = RightModule_Ring_AbelianGroup_Linear_Compatible
    { rightModule_ring           :: RingE
    , rightModule_abelianGroup   :: AbelianGroupE
    , rightModule_actorLinearity :: ActorLinearE
    , rightModule_acteeLinearity :: ActeeLinearE
    , rightModule_compatibility  :: CompatibleE}
  | RightModule_Named Symbol
                     RightModuleE

-- | A right module (v, a) over the ring (r, p, m).
class RightModuleK op p m r a v s

type family RightModuleS (op :: k0) (p :: k1) (m :: k2) r (a :: k3) v :: RightModuleE

type RightModuleC op p m r a v = RightModuleK op p m r a v (RightModuleS op p m r a v)

instance ( RingK p m r zr
         , AbelianGroupK a v zag
         , ActorLinearK 'R m p r a v zor
         , ActeeLinearK 'R m r a v zee
         , CompatibleK 'R op m r v zrc
         ) => RightModuleK op p m r a v
           (RightModule_Ring_AbelianGroup_Linear_Compatible zr zag zor zee zrc)

type RightModule op p m r a v =
  ( RightModuleC op p m r a v
  , Ring p m r
  , AbelianGroup a v
  , LinearActsOn 'R m p r a v
  , RightCompatible op m r v
  )

max_stars_repo_name: evertedsphere/noether

max_issues_repo_name: evertedsphere/noether

max_forks_repo_name: evertedsphere/noether


text:
{-# LANGUAGE ViewPatterns, NoMonomorphismRestriction, FlexibleInstances #-}
module Math.Probably.LinearPDF where

import qualified Math.Probably.Student as S
import Numeric.LinearAlgebra

type PDF a = a->Double

instance Show (a->Double) where
    show f = error "Math.Prob.PDF: showing function"

instance Eq (a->Double) where
    f == g = error "Math.Prob.PDF: comparing functions!"

instance Num a => Num (a->Double) where
    f + g = \x-> f x + g x
    f * g = \x-> f x * g x
    f - g = \x-> f x - g x
    abs f = \x -> abs (f x)
    signum f = \x-> signum $ f x
    fromInteger = error $ "Math.Prob.PDF: frominteger on function"

uniform :: (Real a, Fractional a) => a-> a-> PDF a
uniform from to = \x-> if x>=from && x <=to
                               then realToFrac $ 1/(to-from )
                               else 0

--http://en.wikipedia.org/wiki/Normal_distribution
gauss :: (Real a, Floating a) => a-> a-> PDF a
gauss mean sd = \x->realToFrac $ recip(sd*sqrt(2*pi))*exp(-(x-mean)**2/(2*sd*sd))

gaussD :: Double -> Double-> PDF Double
gaussD mean sd = \x->recip(sd*sqrt(2*pi))*exp(-(x-mean)**2/(2*sd*sd))

gammafun = exp . S.gammaln

--http://en.wikipedia.org/wiki/Gamma_distribution
gammaD :: Double -> Double -> PDF Double
gammaD k theta x = x**(k-1)*(exp(-x/theta)/(theta**k*(exp (S.gammaln k))))

invGammaD :: Double -> Double -> PDF Double
invGammaD a b x =(b**a/gammafun a)*(1/x)**(a+1)*exp(-b/x)

logNormal m sd x = recip (x*sd*sqrt(2*pi)) * exp (negate $ square (log x - m) / (2*square sd))
    where square x = x*x

beta a b x = (recip $ S.beta a b) * x ** (a-1) * (1-x) ** (b-1)
{-# SPECIALIZE gauss :: Double -> Double-> PDF Double #-}

--this is a Prob. Mass function, not density.
binomial :: Int -> Double -> PDF Int
binomial n p k = realToFrac (choose  (toInteger n) (toInteger k)) * p^^k * (1-p)^^(n-k)

zipPdfs :: Num a => PDF a -> PDF b -> PDF (a,b)
zipPdfs fx fy (x,y) = fx x * fy y

mulPdf :: Num a => PDF a -> PDF a -> PDF a
mulPdf d1 d2 = \x -> (d1 x * d2 x)

--instance Num a => Num (PDF a) where

--http://blog.plover.com/math/choose.html
choose :: Integer -> Integer -> Integer
choose n 0 = 1
choose 0 k = 0
choose (n) (k) |  k > n `div` 2 = choose n (n-k)
               |  otherwise = (choose (n-1) (k-1)) * (n) `div` (k)

multiNormal :: Vector Double -> Matrix Double -> PDF (Vector Double)
multiNormal mu sigma =
  let k = realToFrac $ dim mu
      invSigma = inv sigma
      mat1 = head . head . toLists
  in \x-> recip ((2*pi)**(k/2) * sqrt(det sigma)) * exp (mat1 $ negate $ 0.5*(asRow $ x-mu) `multiply` invSigma `multiply` (asColumn $ x-mu) )

mu1 = 2 |> [0, 0::Double]
sig1 = (2><2)[1::Double, 0,
              0, 1]

tst = multiNormal mu1 sig1 mu1
tsta = inv sig1
tstb = det sig1
tstc = let mu = mu1
           sigma = sig1
           x = mu1
           invSigma = inv sigma
       in (asRow $ x-mu) `multiply` invSigma `multiply` (asColumn $ x-mu)

max_stars_repo_name: glutamate/probably

max_issues_repo_name: glutamate/probably

max_forks_repo_name: glutamate/probably


text:
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{-|
Module      : RatioExemplars
Description : Instance for the ContextualizedConcept class for exemplars on a rational measurement scale
Copyright   : (c) Juergen Hahn, 2016
License     : GPL-3
Maintainer  : hahn@geoinfo.tuwien.ac.at
Stability   : tested
-}
module ExemplarScales.RatioExemplars where

import Data.Function (on)
import qualified Data.Vector.Unboxed as U
import Statistics.Sample.KernelDensity (kde)
import qualified Data.List as List
import ContextualizedConcept
import Concept

instance ContextualizedConcept c Double where
  calculatePrototype =calculatePrototype' . createKDE . extractData

calculatePrototype' :: [(Double,Double)] -> (Double, Double)
calculatePrototype' = List.maximumBy (compare `on` snd)

createKDE :: [Double] -> [(Double,Double)]
createKDE  rawdata =U.toList . uncurry U.zip . kde 64 $ dataVector
  where dataVector= U.fromList rawdata

extractData :: Concept c Double -> [Double]
extractData = map getExemplar . toObservationList

max_stars_repo_name: juergenhah/ContextAlgebra

max_issues_repo_name: juergenhah/ContextAlgebra

max_forks_repo_name: juergenhah/ContextAlgebra


text:
module BondCalculus.ODEExtraction
  (IVP(..), PrintStyle(..), ODE(..), AsSage(..),
  matlabExpr, sympyExpr, sageExpr, matlabODE, vectorFieldToODEs, sageODE, sageExprAbst,
  sageExprAbst', extractIVP, sympyODE, solveODEPython, printODEPython,
  sympySimplify, runPython, generateSage) where

import BondCalculus.Base
import BondCalculus.Symbolic
import BondCalculus.Processes
import BondCalculus.Vector
import BondCalculus.AST (pretty)
import qualified BondCalculus.AST as AST
import qualified Data.HashMap.Strict as H
import BondCalculus.Simulation (Trace)
import Data.String.Utils
import Text.Printf

-- import qualified Control.Exception as X
import qualified Data.Map as M
import qualified Data.List as L
import Data.Maybe

--import qualified Numeric.LinearAlgebra as LA

import System.IO.Unsafe
-- import qualified System.Process as OS
import System.Process (readProcess)
-- import Data.Maybe

newtype IVP a = IVP ([String], [Expr (Atom a)], [a]) deriving (Show, Eq)
newtype ODE a = ODE ([String], [Expr (Atom a)])

-------------------------------
-- MATLAB script output:
-------------------------------

extractODE :: (ExprConstant a) => AST.Env -> ConcreteAffinityNetwork a -> P' a -> ODE a
extractODE env network p = vectorFieldToODEs p'
  where -- p' must be expressed in the same basis of p
        -- otherwise the variable order becomes
        -- inconsistient, and we get extra variables
        -- p' :: Vect AST.Species (SymbolicExpr a)
        p' = fromList [(p'coeff i, i) | (_,i) <- toList p]
        -- p'coeff :: AST.Species -> SymbolicExpr a
        p'coeff i = fromMaybe 0.0 (H.lookup i v)
        Vect v = dPdt' tr network p
        tr = tracesGivenNetwork network env

extractIVP :: ExprConstant a => AST.Env -> ConcreteAffinityNetwork a -> P' a -> [a] -> IVP a
extractIVP env network p = fromODEToIVP (extractODE env network p)

vectorFieldToODEs :: ExprConstant a => P' a -> ODE a
vectorFieldToODEs v = ODE (vars, rhss)
  where
    vars = map (pretty.snd) kis
    rhss = map (simplify.fst) kis
    kis = toList v

fromODEToIVP :: ODE a -> [a] -> IVP a
fromODEToIVP (ODE (x,y)) z = IVP (x, y, take (length x) z)

-- modelToIVP :: BondCalculusModel SymbolicExpr Double -> IVP

matlabODE :: IVP Double -> (Int,(Double,Double)) -> Either String String
matlabODE ivp (n,(t0,tn)) =
  if any isNothing eqns
  then Left "An ODE equation has an unbound variable"
  else Right $
  "init = [" ++ L.intercalate ";" (map show inits) ++ "];\n" ++
    "t = linspace(" ++ show t0 ++ "," ++ show tn ++ "," ++ show n ++ ");\n" ++
    "function xdot = f(x,t) \n\n" ++
    L.intercalate ";\n" (catMaybes eqns) ++
    ";\nendfunction;\n\n" ++
    -- Need new support for calculating jacobian
    -- "function jac = jj(x,t)\n\n" ++
    -- matlabJac env p' ++ ";\nendfunction;\n\n" ++
    "x = lsode(f, init, t);\n" ++
    "save (\"-ascii\", \"-\", \"x\");"
      where
        xdots = ["xdot(" ++ show i ++ ")" | i <- [(1::Integer)..]]
        xvars = ["x(" ++ show i ++ ")" | i <- [(1::Integer)..]]
        varmap = M.fromList $ zip vars xvars
        IVP (vars,rhss,inits) = ivp
        eqns = [fmap (\y -> xdot ++ " = " ++ y) (matlabExpr varmap rhs)
               | (xdot,rhs) <- zip xdots rhss]

class Show a => AsSage a where
    asSage :: a -> String

instance AsSage Double where
    asSage = show

instance AsSage Interval where
    asSage x = printf "RIF('[%.20f .. %.20f]')" l u
        where (l, u) = endpoints x
              -- Do some fiddling with the endpoints to ensure sound enclosure
            --   l' = l - 6e-21
            --   u' = u + 6e-21

sageExpr :: AsSage a => M.Map String String -> Expr (Atom a) -> Maybe String
sageExpr mp (Atom (Var x)) = M.lookup x mp
sageExpr mp (Atom (Const x)) = return $ asSage x
sageExpr mp (x `Sum` y) = do
    x' <- sageExpr mp x
    y' <- sageExpr mp y
    return $ x' ++ " + " ++ y'
sageExpr mp (x `Prod` y) = do
    x' <- sageExpr mp x
    y' <- sageExpr mp y
    return $ "(" ++ x' ++ ") * (" ++ y' ++ ")"
sageExpr mp (x `Frac` y) = do
    x' <- sageExpr mp x
    y' <- sageExpr mp y
    return $ "(" ++ x' ++ ") / (" ++ y' ++ ")"
sageExpr mp (x `Pow` y) = do
    x' <- sageExpr mp x
    y' <- sageExpr mp y
    return $ "(" ++ x' ++ ") ** (" ++ y' ++ ")"
sageExpr mp (Abs x) = do
    x' <- sageExpr mp x
    return $ "abs(" ++ x' ++ ")"
sageExpr mp (Sin x) = do
    x' <- sageExpr mp x
    return $ "sin(" ++ x' ++ ")"
sageExpr mp (Cos x) = do
    x' <- sageExpr mp x
    return $ "sin(" ++ x' ++ ")"
sageExpr mp (Tan x) = do
    x' <- sageExpr mp x
    return $ "tan(" ++ x' ++ ")"
sageExpr mp (Exp x) = do
    x' <- sageExpr mp x
    return $ "exp(" ++ x' ++ ")"
sageExpr mp (Log x) = do
    x' <- sageExpr mp x
    return $ "log(" ++ x' ++ ")"

-- Convert a sage expression into a string with all non-simple double constants (e.g. intervals)
-- abstracted by symbolic variables
sageExprAbst' :: (AsSage a, Boundable a) => M.Map String String -> Expr (Atom a) -> Maybe (M.Map String String, String)
sageExprAbst' = sageExprAbst M.empty

sageExprAbst :: (AsSage a, Boundable a) => M.Map String String -> M.Map String String -> Expr (Atom a) -> Maybe (M.Map String String, String)
sageExprAbst cmp mp (Atom (Var x)) = do
    x' <- M.lookup x mp
    return (cmp, x')
sageExprAbst cmp mp (Atom (Const x)) = return $ case singleValue x of
    Just 1.0    -> (cmp, "1")
    Just (-1.0) -> (cmp, "-1")
    Just 0.0    -> (cmp, "0")
    Just v      -> (cmp, asSage v)
    Nothing     -> (cmp', s)
        where n = M.size cmp
              s = "a" ++ show n
              cmp' = M.insert s (asSage x) cmp
sageExprAbst cmp mp (x `Sum` y) = f2 cmp mp x y $ \x' y' -> x' ++ " + " ++ y'
sageExprAbst cmp mp (x `Prod` y) = f2 cmp mp x y $ \x' y' -> "(" ++ x' ++ ") * (" ++ y' ++ ")"
sageExprAbst cmp mp (x `Frac` y) = f2 cmp mp x y $ \x' y' -> "(" ++ x' ++ ") / (" ++ y' ++ ")"
sageExprAbst cmp mp (x `Pow` y) = f2 cmp mp x y $ \x' y' -> "(" ++ x' ++ ") ** (" ++ y' ++ ")"
sageExprAbst cmp mp (Abs x) = f1 cmp mp x $ \x' -> "abs(" ++ x' ++ ")"
sageExprAbst cmp mp (Sin x) = f1 cmp mp x $ \x' -> "sin(" ++ x' ++ ")"
sageExprAbst cmp mp (Cos x) = f1 cmp mp x $ \x' -> "cos(" ++ x' ++ ")"
sageExprAbst cmp mp (Tan x) = f1 cmp mp x $ \x' -> "tan(" ++ x' ++ ")"
sageExprAbst cmp mp (Exp x) = f1 cmp mp x $ \x' -> "exp(" ++ x' ++ ")"
sageExprAbst cmp mp (Log x) = f1 cmp mp x $ \x' -> "log(" ++ x' ++ ")"
sageExprAbst cmp mp x = error $ "Support for converting expression " ++ show x ++ " to sage/sympy not implemented!"

f2 cmp mp x y g = do
    (cmp', x') <- sageExprAbst cmp mp x
    (cmp'', y') <- sageExprAbst cmp' mp y
    return (cmp'', g x' y')

f1 cmp mp x g = do
    (cmp', x') <- sageExprAbst cmp mp x
    return (cmp', g x')

-- This assumes we can reduce the ODEs to sage
sageODE :: (AsSage a, Boundable a) => IVP a -> String -> Either String String
sageODE (IVP (vars, rhss, inits)) network = case odeExprs of
        Nothing -> Left "An ODE equation has an unbound variable"
        Just (cmp, odes) -> Right $
            "from sage.all import *\n" ++
            "import sympy as sym\n\n" ++
            "R, x = PolynomialRing(RIF, " ++ show nvars ++ ", ', '.join(map('x{}'.format, range(0, " ++ show nvars ++ ")))).objgens()\n" ++
            "xsr = [SR.var(str(x1)) for x1 in x]\n" ++
            "for v in xsr:\n"++
            "   assume(v, 'real')\n"++
            "xsym = sym.var(','.join(map('x{}'.format, range(0," ++ show nvars ++ "))))\n" ++
            (if (M.size cmp) > 0
             then "asym = sym.var(','.join(map('a{}'.format, range(0," ++ show (M.size cmp) ++ "))))\n"
             else "asym = []\n") ++
            "y0 = [" ++ L.intercalate ", " (map asSage inits) ++ "]\n" ++
            "ysymraw = [" ++ L.intercalate ", " odes ++ "]\n" ++
            "ysym = [sym.simplify(y1) for y1 in ysymraw]\n" ++
            "ysr = [y1._sage_().substitute(" ++ subsExpr ++ ") for y1 in ysym]\n" ++
            "varmap = {" ++ L.intercalate ", " varmappings ++ "}\n" ++
            "varmapsr = {" ++ L.intercalate ", " varmappingssr ++ "}\n" ++
            "affinity_network = \"" ++ network ++ "\"\n" ++
            "try:\n" ++
            "    y = vector([R(y1) for y1 in ysr])\n" ++
            "except TypeError:\n" ++
            "    y = None\n"
            where subsExpr = L.intercalate ","
                             [x ++ "=" ++ y | (x, y) <- M.toList cmp]
                  varmappings = [ "'" ++ p ++ "': x[" ++ show i ++ "]"
                                | p <- vars | i <- [0..] ]
                  varmappingssr = [ "'" ++ p ++ "': xsr[" ++ show i ++ "]"
                                  | p <- vars | i <- [0..] ]
    where
    --   odeExprs' :: AsSage a => M.Map String String -> [Expr (Atom a)] -> Maybe (M.Map String String, [String])
      odeExprs' cmp (x:xs) = do
        (cmp', expr) <- sageExprAbst cmp varmap x
        (cmp'', exprs) <- odeExprs' cmp' xs
        return (cmp'', expr:exprs)
      odeExprs' cmp [] = Just (cmp, [])
      odeExprs :: Maybe (M.Map String String, [String])
      odeExprs = odeExprs' M.empty rhss
      nvars = length vars
      xvars = ["x" ++ show n | (n, _) <- zip [0..] vars]
      varmap = M.fromList $ zip vars xvars
    --   odes = map odeExpr rhss
    --   odeExpr rhs = do
    --     expr <- sageExpr varmap rhs
    --     return $ "sym.simplify(" ++ expr ++ ")"

sympyODE :: IVP Double -> (Int,(Double,Double)) -> Either String String
sympyODE ivp (n,(t0,tn)) =
  if any isNothing eqns
  then Left "An ODE equation has an unbound variable"
  else Right $
  "import sympy as sym\n" ++
  "import numpy as np\n" ++
  "from sympy.abc import t\n" ++
  "from scipy.integrate import odeint\n" ++
  "import sys\n" ++
  "sys.setrecursionlimit(100000)\n\n" ++
  "xs = [" ++ L.intercalate ","
    ["sym.Function('" ++ show x ++ "')" | x <- vars ] ++ "]\n" ++
  "xts = [x(t) for x in xs]\n" ++
  "odes = [" ++ L.intercalate ", " (catMaybes eqns) ++ "]\n" ++
  "y0 = [" ++ L.intercalate ", " (map show inits) ++ "]\n" ++
  "ts = np.linspace(" ++ show t0 ++ "," ++ show tn ++ "," ++ show n ++ ")\n" ++
  "rhss = [eqn.rhs for eqn in odes]\n" ++
  "Jac = sym.Matrix(rhss).jacobian(xts)\n" ++
  "f = sym.lambdify((xts, t), rhss, modules='numpy')\n" ++
  "J = sym.lambdify((xts, t), Jac, modules='numpy')\n" ++
  -- Try to solve ODEs using symbolic Jacobian
  "try: ys = odeint(f, y0, ts, (), J)\n" ++
  -- Fallback to using method without Jacobian
  "except NameError: ys = odeint(f, y0, ts, ())\n" ++
  -- "print(ys)\n" ++
  "print('\\n'.join([' '.join([('%.18e' % y).replace('nan', '0') for y in ysa]) for ysa in ys]))"
    where
      xvars = ["xts[" ++ show i ++ "]" | i <- [(0::Integer)..]]
      xdots = [xvar ++ ".diff()" | xvar <- xvars]
      varmap = M.fromList $ zip vars xvars
      IVP (vars,rhss,inits) = ivp
      mkeqn xdot y = "sym.Eq(" ++ xdot ++ ", sym.simplify(" ++ y ++ "))"
      eqns = [fmap (mkeqn xdot) (sympyExpr varmap rhs)
             | (xdot,rhs) <- zip xdots rhss]

data PrintStyle = Plain | Pretty | LaTeX | MathML deriving (Eq)

sympyODEPrint :: ODE Double -> PrintStyle -> Either String String
sympyODEPrint ode style =
  if any isNothing eqns
  then Left "An ODE equation has an unbound variable"
  else Right $
  "import sympy as sym\n" ++
  "import numpy as np\n" ++
  "from sympy.abc import t\n" ++
  "import sys\n" ++
  "sys.setrecursionlimit(100000)\n\n" ++
  -- "sym.init_printing(" ++ printingoptions ++ ")\n\n" ++
  "names = [" ++ L.intercalate "," ["'[" ++ (if style == LaTeX then replace "->" "\\\\rightarrow " (replace " " "\\\\ " x) else x) ++ "]'" | x <- vars] ++ "]\n" ++
  "xs = list(map(sym.Function, names))\n" ++
  "xts = [x(t) for x in xs]\n" ++
  "odes = [" ++ L.intercalate ", " (catMaybes eqns) ++ "]\n" ++
  "rhss = [eqn.rhs for eqn in odes]\n" ++
  case style of
    Pretty -> "for ode in odes: print(sym.pretty(ode))\n"
    Plain -> "for ode in odes: print(ode)\n"
    -- print system of odes as latex (formatted for use in an align environment)
    LaTeX -> "print(sym.latex(odes)[7:-7].replace('=', '& =').replace(', \\\\  ', ' \\\\\\\\\\n').replace('{\\\\left(t \\\\right)}', '').replace('\\\\frac{d}{d t}', '\\\\frac{\\\\mathrm d}{\\\\mathrm d t}').replace('operatorname', 'mathrm'))\n"
    MathML -> "from sympy.printing import print_mathml\nprint_mathml(odes)\n"
    where
      xvars = ["xts[" ++ show i ++ "]" | i <- [(0::Integer)..]]
      xdots = [xvar ++ ".diff()" | xvar <- xvars]
      varmap = M.fromList $ zip vars xvars
      ODE (vars,rhss) = ode
      mkeqn xdot y = "sym.Eq(" ++ xdot ++ ", sym.simplify(" ++ y ++ "))"
      eqns = [fmap (mkeqn xdot) (sympyExpr varmap rhs)
             | (xdot,rhs) <- zip xdots rhss]

matlabExpr :: M.Map String String -> SymbolicExpr Double -> Maybe String
matlabExpr mp (Atom (Var x)) = M.lookup x mp
matlabExpr _ (Atom (Const x)) = return $ show x
matlabExpr mp (x `Sum` y) = do
  x' <- matlabExpr mp x
  y' <- matlabExpr mp y
  return $ x' ++ " .+ " ++ y'
matlabExpr mp (x `Prod` y) = do
  x' <- matlabExpr mp x
  y' <- matlabExpr mp y
  return $ "(" ++ x' ++ ") .* (" ++ y' ++ ")"
matlabExpr mp (x `Pow` y) = do
  x' <- matlabExpr mp x
  y' <- matlabExpr mp y
  return $ "(" ++ x' ++ ") .** (" ++ y' ++ ")"
matlabExpr mp (Abs x) = do
  x' <- matlabExpr mp x
  return $ "abs(" ++ x' ++ ")"

sympyExpr :: M.Map String String -> SymbolicExpr Double -> Maybe String
sympyExpr mp (Atom (Var x)) = M.lookup x mp
sympyExpr _ (Atom (Const x)) = return $ show x
sympyExpr mp (x `Sum` y) = do
  x' <- sympyExpr mp x
  y' <- sympyExpr mp y
  return $ x' ++ " + " ++ y'
sympyExpr mp (x `Prod` y) = do
  x' <- sympyExpr mp x
  y' <- sympyExpr mp y
  return $ "(" ++ x' ++ ") * (" ++ y' ++ ")"
sympyExpr mp (x `Frac` y) = do
  x' <- sympyExpr mp x
  y' <- sympyExpr mp y
  return $ "(" ++ x' ++ ") / (" ++ y' ++ ")"
sympyExpr mp (x `Pow` y) = do
  x' <- sympyExpr mp x
  y' <- sympyExpr mp y
  return $ "(" ++ x' ++ ") ** (" ++ y' ++ ")"
sympyExpr mp (Abs x) = do
  x' <- sympyExpr mp x
  -- return $ x'
  return $ "abs(" ++ x' ++ ")"
sympyExpr mp (Sin x) = do
  x' <- sympyExpr mp x
  return $ "sin(" ++ x' ++ ")"
sympyExpr mp (Cos x) = do
  x' <- sympyExpr mp x
  return $ "sin(" ++ x' ++ ")"
sympyExpr mp (Tan x) = do
  x' <- sympyExpr mp x
  return $ "tan(" ++ x' ++ ")"
sympyExpr mp (Exp x) = do
  x' <- sympyExpr mp x
  return $ "exp(" ++ x' ++ ")"
sympyExpr mp (Log x) = do
  x' <- sympyExpr mp x
  return $ "log(" ++ x' ++ ")"

---------------------------------------
-- Using Octave to execute the scripts
---------------------------------------

runPython :: String -> IO String
runPython = readProcess "python3" ["-q"]

callSolveODEPython :: AST.Env -> ConcreteAffinityNetwork Double -> P' Double -> [Double] -> (Int, (Double, Double)) -> IO String
callSolveODEPython env network p inits tr = case sympyODE (extractIVP env network p inits) tr of
  Right script -> do
    putStrLn $ "Python script:\n\n" ++ script
    writeFile "script.py" script
    runPython script
  Left _ -> undefined

generateSage :: (ExprConstant a, AsSage a)
             => String
             -> AST.Env
             -> ConcreteAffinityNetwork a
             -> P' a
             -> [a]
             -> IO String
generateSage filename env network p inits = case sageODE (extractIVP env network p inits) (pretty network) of
  Right script -> do
    putStrLn $ "Sage script:\n\n" ++ script
    writeFile filename script
    return script
  Left _ -> undefined

callPrintODEPython :: AST.Env -> ConcreteAffinityNetwork Double -> P' Double -> PrintStyle -> IO String
callPrintODEPython env network p style = case sympyODEPrint (extractODE env network p) style of
  Right script -> do
    -- putStrLn $ "Python script:\n\n" ++ script
    writeFile "print_script.py" script
    runPython script
  Left _ -> undefined

solveODEPython :: AST.Env -> ConcreteAffinityNetwork Double -> P' Double -> [Double] -> (Int, (Double, Double)) -> Trace
solveODEPython env network p inits tr@(n,(t0,tn))
  = let raw = unsafePerformIO (callSolveODEPython env network p inits tr)
        ts = [t0 + fromIntegral i*(tn-t0)/fromIntegral n | i <- [0..n]]
        yss = (map (map read.words) $ lines raw) :: [[Double]]
        ys = [fromList (xs `zip` pbasis) | xs <- yss]
        pbasis = map snd $ toList p
    in ts `zip` ys

sympySimplify :: String -> String
sympySimplify s = unsafePerformIO (runPython script)
  where script = "from sympy import simplify, sympify\nimport re\n" ++
                 "print(simplify(sympify(re.sub(r'([0-9]+)\\.0(?![0-9]*[1-9])', r'\\1', " ++ show s ++ "))))"

printODEPython :: AST.Env -> ConcreteAffinityNetwork Double -> P' Double -> PrintStyle -> String
printODEPython env network p style
  = let raw = unsafePerformIO (callPrintODEPython env network p style)
    in raw
-- callOctave env p mts p' ts = let
--     script = matlabODE env (wholeProc env p mts) p' ts
--   in do
--     putStrLn $ "Octave Script: \n" ++ script
--     OS.readProcess
--       "octave" ["-q", "--eval", script] []

-- | Solver which calculates the symbolic Jacobian, writes MATLAB code, and executes it with GNU Octave. (General purpose, deals with stiff systems, uses LSODE.)
-- solveODEoctave :: Solver
-- solveODEoctave env p mts p' ts@(n,(t0,tn))
--     = let raw = unsafePerformIO (callOctave env p mts p' ts)
--       in (n>< Map.size p') $ map s2d $ words raw

-- Return the MATLAB script for ODEs
-- matlabScript :: Env
--              -> Process
--              -> MTS
--              -> P'
--              -> (Int, (Double, Double))
--              -> String
-- matlabScript env p mts = matlabODE env (wholeProc env p mts)

max_stars_repo_name: twright/bondwb

max_issues_repo_name: twright/bondwb

max_forks_repo_name: twright/bondwb


text:
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE FlexibleInstances   #-}
{-# LANGUAGE PolyKinds           #-}
{-# LANGUAGE ScopedTypeVariables #-}
-- |
-- Module      : Test.ShowType
-- Copyright   : [2017..2020] The Accelerate Team
-- License     : BSD3
--
-- Maintainer  : Trevor L. McDonell <trevor.mcdonell@gmail.com>
-- Stability   : experimental
-- Portability : non-portable (GHC extensions)
--

module Test.ShowType
  where

import Data.Complex

data ArgType (a :: *) = AT

showType :: forall proxy a. Show (ArgType a) => proxy a -> String
showType _ = show (AT :: ArgType a)

instance Show (ArgType a) => Show (ArgType (Complex a)) where
  show _ = "Complex " ++ show (AT :: ArgType a)

instance Show (ArgType Float)  where show _ = "Float"
instance Show (ArgType Double) where show _ = "Double"

max_stars_repo_name: tmcdonell/accelerate-fft

max_issues_repo_name: tmcdonell/accelerate-fft

max_forks_repo_name: tmcdonell/accelerate-fft