text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Kalman
-- Copyright   :  (c) 2016 FP Complete Corporation
-- License     :  MIT (see LICENSE)
-- Maintainer  :  dominic@steinitz.org
--
-- =The Theory
-- The model for the extended Kalman filter is given by
--
-- \[
-- \begin{aligned}
-- \boldsymbol{x}_i &= \boldsymbol{a}_i(\boldsymbol{x}_{i-1}, \boldsymbol{u}_{i-1}) + \boldsymbol{\psi}_{i-1} \\
-- \boldsymbol{y}_i &= \boldsymbol{h}_i(\boldsymbol{x}_i) + \boldsymbol{\upsilon}_i
-- \end{aligned}
-- \]
--
-- where
--
-- * \(\boldsymbol{a_i}\)\ is some non-linear vector-valued possibly
-- time-varying state update function.
--
-- * \(\boldsymbol{x_i}\) is the state at time \(i\).
--
-- * \(\boldsymbol{u_i}\) is the control at time \(i\).
--
-- * \(\boldsymbol{\psi}_{i}\) are independent normally
-- distributed random variables with mean 0 representing the fact that
-- the state update is noisy: \(\boldsymbol{\psi}_{i} \sim {\cal{N}}(0,Q_i)\).
--
-- * \(\boldsymbol{h}_i\)\ is some non-linear vector-valued possibly time-varying
-- function describing how we observe the hidden state process.
--
-- * \(\boldsymbol{y_i}\) is the observable at time \(i\).
--
-- * \(\boldsymbol{\upsilon}_i\) are independent normally
-- distributed random variables with mean 0 represent the fact that
-- the observations are noisy: \(\boldsymbol{\upsilon}_{i} \sim {\cal{N}}(0,R_i)\).
--
-- Note that in most presentations of the Kalman filter (the
-- [wikipedia](https://en.wikipedia.org/wiki/Kalman_filter)
-- presentation being an exception), the state update function and the
-- observation function are taken to be constant over time as is the
-- noise for the state and the noise for
-- the observation. In symbols, \(\forall i \, \boldsymbol{a}_i = \boldsymbol{a}\)
-- for some \(\boldsymbol{a}\), \(\forall i \, \boldsymbol{h}_i = \boldsymbol{h}\)
-- for some \(\boldsymbol{h}\),
-- \(\forall i \, Q_i = Q\) for some \(Q\) and \(\forall i \, \boldsymbol{a}_i = R\)
-- for some \(R\).
--
-- We assume the whole process starts at 0 with our belief of the state (aka the
-- prior state) being given by
-- \(\boldsymbol{x}_0 \sim {\cal{N}}(\boldsymbol{\mu}_0, \Sigma_0)\)
--
-- The Kalman filtering process is a recursive procedure as follows:
--
-- (1) Make a prediction of the current system state, given our previous
-- estimation of the system state.
--
-- (2) Update our prediction, given a newly acquired measurement.
--
-- The prediction and update steps depend on our system and measurement
-- models, as well as our current estimate of the system state.
--
-- = An Extended Example
--
-- The equation of motion for a pendulum of unit length subject to
-- [Gaussian white
-- noise](https://en.wikipedia.org/wiki/White_noise#Mathematical_definitions)
-- is
--
-- \[
-- \frac{\mathrm{d}^2\alpha}{\mathrm{d}t^2} = -g\sin\alpha + w(t)
-- \]
--
-- We can discretize this via the usual [Euler method](https://en.wikipedia.org/wiki/Euler_method)
--
-- \[
-- \begin{bmatrix}
-- x_{1,i} \\
-- x_{2,i}
-- \end{bmatrix}
-- =
-- \begin{bmatrix}
-- x_{1,i-1} + x_{2,i-1}\Delta t \\
-- x_{2,i-1} - g\sin x_{1,i-1}\Delta t
-- \end{bmatrix}
-- +
-- \mathbf{q}_{i-1}
-- \]
--
-- where \(q_i \sim {\mathcal{N}}(0,Q)\) and
--
-- \[
-- Q
-- =
-- \begin{bmatrix}
-- \frac{q^c \Delta t^3}{3} & \frac{q^c \Delta t^2}{2} \\
-- \frac{q^c \Delta t^2}{2} & {q^c \Delta t}
-- \end{bmatrix}
-- \]
--
-- Assume that we can only measure the horizontal position of the
-- pendulum and further that this measurement is subject to error so that
--
-- \[
-- y_i = \sin x_i + r_i
-- \]
--
-- where \(r_i \sim {\mathcal{N}}(0,R)\).
--
-- First let's set the time step and the acceleration caused by earth's gravity.
--
-- > {-# LANGUAGE DataKinds #-}
-- >
-- > import Numeric.Kalman
-- >
-- > deltaT, g :: Double
-- > deltaT = 0.01
-- > g  = 9.81
-- > bigQ :: Sym 2
-- > bigQ = sym $ matrix bigQl
-- >
-- > qc :: Double
-- > qc = 0.01
-- >
-- > bigQl :: [Double]
-- > bigQl = [ qc * deltaT^3 / 3, qc * deltaT^2 / 2,
-- >           qc * deltaT^2 / 2,       qc * deltaT
-- >         ]
-- >
-- > bigR :: Sym 1
-- > bigR  = sym $ matrix [0.1]
-- > stateUpdate :: R 2 -> R 2
-- > stateUpdate u =  vector [x1 + x2 * deltaT, x2 - g * (sin x1) * deltaT]
-- >   where
-- >     (x1, w) = headTail u
-- >     (x2, _) = headTail w
-- >
-- > observe :: R 2 -> R 1
-- > observe a = vector [sin x] where x = fst $ headTail a
-- > linearizedObserve :: R 2 -> L 1 2
-- > linearizedObserve a = matrix [cos x, 0.0] where x = fst $ headTail a
-- >
-- > linearizedStateUpdate :: R 2 -> Sq 2
-- > linearizedStateUpdate u = matrix [1.0,                    deltaT,
-- >                                   -g * (cos x1) * deltaT,    1.0]
-- >   where
-- >     (x1, _) = headTail u
--
-- Now we can create extended and unscented filters which consume a
-- single observation.
--
-- > singleEKF = runEKF (const observe) (const linearizedObserve) (const bigR)
-- >              (const stateUpdate) (const linearizedStateUpdate) (const bigQ)
-- >              undefined
--
-- > singleUKF = runUKF (const observe) (const bigR) (const stateUpdate) (const bigQ)
-- >              undefined
--
-- We start off not too far from the actual value.
--
-- > initialDist =  (vector [1.6, 0.0],
-- >                 sym $ matrix [0.1, 0.0,
-- >                               0.0, 0.1])
--
-- Using some data generated using code made available with Simo
-- Särkkä's
-- [book](http://www.cambridge.org/gb/academic/subjects/statistics-probability/applied-probability-and-stochastic-networks/bayesian-filtering-and-smoothing),
-- we can track the pendulum using the extended Kalman filter.
--
-- > multiEKF obs = scanl singleEKF initialDist (map (vector . pure) obs)
--
-- And then plot the results.
--
-- <<diagrams/src_Numeric_Kalman_diagE.svg#diagram=diagE&height=600&width=500>>
--
-- And also track it using the unscented Kalman filter.
--
-- > multiUKF obs = scanl singleUKF initialDist (map (vector . pure) obs)
--
-- And also plot the results
--
-- <<diagrams/src_Numeric_Kalman_diagU.svg#diagram=diagU&height=600&width=500>>
--
-- === Code for Plotting
--
-- The full code for plotting the results:
--
-- > import qualified Graphics.Rendering.Chart as C
-- > import Graphics.Rendering.Chart.Backend.Diagrams
-- > import Data.Colour
-- > import Data.Colour.Names
-- > import Data.Default.Class
-- > import Control.Lens
-- >
-- > import Data.Csv
-- > import System.IO hiding ( hGetContents )
-- > import Data.ByteString.Lazy ( hGetContents )
-- > import qualified Data.Vector as V
-- >
-- > import Numeric.LinearAlgebra.Static
-- >
-- > chartEstimated :: String ->
-- >               [(Double, Double)] ->
-- >               [(Double, Double)] ->
-- >               [(Double, Double)] ->
-- >               C.Renderable ()
-- > chartEstimated title acts obs ests = C.toRenderable layout
-- >   where
-- >
-- >     actuals = C.plot_lines_values .~ [acts]
-- >             $ C.plot_lines_style  . C.line_color .~ opaque red
-- >             $ C.plot_lines_title .~ "Actual Trajectory"
-- >             $ C.plot_lines_style  . C.line_width .~ 1.0
-- >             $ def
-- >
-- >     measurements = C.plot_points_values .~ obs
-- >                  $ C.plot_points_style  . C.point_color .~ opaque blue
-- >                  $ C.plot_points_title .~ "Measurements"
-- >                  $ def
-- >
-- >     estimas = C.plot_lines_values .~ [ests]
-- >             $ C.plot_lines_style  . C.line_color .~ opaque black
-- >             $ C.plot_lines_title .~ "Inferred Trajectory"
-- >             $ C.plot_lines_style  . C.line_width .~ 1.0
-- >             $ def
-- >
-- >     layout = C.layout_title .~ title
-- >            $ C.layout_plots .~ [C.toPlot actuals, C.toPlot measurements, C.toPlot estimas]
-- >            $ C.layout_y_axis . C.laxis_title .~ "Angle / Horizontal Displacement"
-- >            $ C.layout_y_axis . C.laxis_override .~ C.axisGridHide
-- >            $ C.layout_x_axis . C.laxis_title .~ "Time"
-- >            $ C.layout_x_axis . C.laxis_override .~ C.axisGridHide
-- >            $ def
-- >
-- > diagE = do
-- >   h <- openFile "data/matlabRNGs.csv" ReadMode
-- >   cs <- hGetContents h
-- >   let df = (decode NoHeader cs) :: Either String (V.Vector (Double, Double))
-- >   case df of
-- >     Left _ -> error "Whatever"
-- >     Right generatedSamples -> do
-- >       let xs = take 500 (multiEKF $ V.toList $ V.map fst generatedSamples)
-- >       let mus = map (fst . headTail . fst) xs
-- >       let obs = V.toList $ V.map fst generatedSamples
-- >       let acts = V.toList $ V.map snd generatedSamples
-- >       denv <- defaultEnv C.vectorAlignmentFns 600 500
-- >       let charte = chartEstimated "Extended Kalman Filter"
-- >                                   (zip [0,1..] acts)
-- >                                   (zip [0,1..] obs)
-- >                                   (zip [0,1..] mus)
-- >       return $ fst $ runBackend denv (C.render charte (600, 500))
-- >
-- > diagU = do
-- >   h <- openFile "data/matlabRNGs.csv" ReadMode
-- >   cs <- hGetContents h
-- >   let df = (decode NoHeader cs) :: Either String (V.Vector (Double, Double))
-- >   case df of
-- >     Left _ -> error "Whatever"
-- >     Right generatedSamples -> do
-- >       let ys = take 500 (multiUKF $ V.toList $ V.map fst generatedSamples)
-- >       let nus = map (fst . headTail . fst) ys
-- >       let obs = V.toList $ V.map fst generatedSamples
-- >       let acts = V.toList $ V.map snd generatedSamples
-- >       denv <- defaultEnv C.vectorAlignmentFns 600 500
-- >       let charte = chartEstimated "Unscented Kalman Filter"
-- >                                   (zip [0,1..] acts)
-- >                                   (zip [0,1..] obs)
-- >                                   (zip [0,1..] nus)
-- >       return $ fst $ runBackend denv (C.render charte (600, 500))
--
-----------------------------------------------------------------------------

module Numeric.Kalman
       (
         runKF, runKFPrediction, runKFUpdate
       , runEKF, runEKFPrediction, runEKFUpdate
       , runUKF, runUKFPrediction, runUKFUpdate
       , runKS, runEKS, runUKS
       )
       where

import GHC.TypeLits
import Numeric.LinearAlgebra.Static
import Data.Maybe ( fromJust )
import Prelude hiding ((<>))

-- | Given our system model and our previous estimate of the system state,
-- we generate a prediction of the current system state by taking the
-- mean of our estimated state as our point estimate and evolving it one
-- step forward with the system evolution function.  In other words, we
-- predict that the current system state is the minimum mean squared
-- error (MMSE) state, which corresponds to the maximum a posteriori
-- (MAP) estimate, from the previous estimate evolved one step by the
-- system evolution function of our system model.
--
-- We also generate a predicted covariance matrix by updating the current
-- covariance and adding any noise estimate.  Updating the current
-- covariance requires the linearized version of the system evolution,
-- consequently the deviation of the actual system evolution from its
-- linearization is manifested in the newly predicted covariance.
--
-- Taken together, the prediction step generates a new Gaussian
-- distribution of predicted system states with mean the evolved MAP
-- state from the previous estimate.
runEKFPrediction
  :: KnownNat n
  => (a -> R n -> R n)  -- ^ System evolution function
                        -- \(\boldsymbol{a}_i(\boldsymbol{u}, \boldsymbol{x})\),
                        -- note the order of the control and the
                        -- input
  -> (a -> R n -> Sq n) -- ^ Linearization of the system evolution at a point \(\frac{\partial \boldsymbol{a}_i}{\partial \boldsymbol{x}}\big|_{{\boldsymbol{u}}, \, {\boldsymbol{x}}}\)
  -> (a -> Sym n)       -- ^ Covariance matrix encoding system evolution noise \(Q_i\)
  -> a                  -- ^ Dynamical input \(\boldsymbol{u}_{i-1}\)
  -> (R n, Sym n)       -- ^ Current estimate \((\hat{\boldsymbol{x}}_{n-1}, \hat{\boldsymbol{\Sigma}}_{n-1})\)
  -> (R n, Sym n)       -- ^ New prediction \(({\boldsymbol{x}_i^\flat}, {\boldsymbol{\Sigma}_i^\flat})\)

runEKFPrediction evolve linEvol sysCov input (estMu, estCov) =
  (predMu, predCov)
  where
    predMu  = evolve input estMu
    predCov = sym $ lin <> unSym estCov <> tr lin + unSym (sysCov input)
    lin     = linEvol input estMu

runKFPrediction
  :: (KnownNat n)
  => (a -> R n -> Sq n) -- ^ Linear system evolution at a point
                        -- \(\boldsymbol{a}_i(\boldsymbol{u}_{i-1}, \boldsymbol{x}_{i-1})\),
                        -- note the order of the control and the
                        -- input
  -> (a -> Sym n)       -- ^ Covariance matrix encoding system evolution noise \(Q_i\)
  -> a                  -- ^ Dynamical input \(\boldsymbol{u}_{i-1}\)
  -> (R n, Sym n)       -- ^ Current estimate
                        -- \((\hat{\boldsymbol{x}}_{i-1}, \hat{\boldsymbol{\Sigma}}_{i-1})\)
  -> (R n, Sym n)       -- ^ New prediction
                        -- \(({\boldsymbol{x}_i^\flat}, {\boldsymbol{\Sigma}_i^\flat})\)
runKFPrediction linEvol =
  runEKFPrediction (\inp sys -> linEvol inp sys #> sys) linEvol

-- | After a new measurement has been taken, we update our original
-- prediction of the current system state using the result of this
-- measurement.  This is accomplished by looking at how much the
-- measurement deviated from our prediction, and then updating our
-- estimated state accordingly.
--
-- This step requires the linearized measurement transformation.
runEKFUpdate
  :: (KnownNat m, KnownNat n)
  => (a -> R n -> R m)   -- ^ System measurement function \(\boldsymbol{h}_i\)
  -> (a -> R n -> L m n) -- ^ Linearization of the observation at a point
                         -- \(\frac{\partial \boldsymbol{h}_i}{\partial \boldsymbol{x}}\big|_{\boldsymbol{x}}\)
  -> (a -> Sym m)        -- ^ Covariance matrix encoding measurement noise \(R_i\)
  -> a                   -- ^ Dynamical input \(\boldsymbol{u}_{i-1}\)
  -> (R n, Sym n)        -- ^ Current prediction \(({\boldsymbol{x}}_i^\flat, {\boldsymbol{\Sigma}}_i^\flat)\)
  -> R m                 -- ^ New measurement \(\boldsymbol{y}_i\)
  -> (R n, Sym n)        -- ^ Updated prediction \((\hat{\boldsymbol{x}}_i, \hat{\boldsymbol{\Sigma}}_i)\)

runEKFUpdate measure linMeas measCov input (predMu, predCov') newMeas =
  (newMu, newCov)
  where
    newMu  = predMu + kkMat #> voff
    newCov = sym $ predCov - kkMat <> skMat <> tr kkMat

    predCov = unSym predCov'

    lin   = linMeas input predMu
    voff  = newMeas - measure input predMu
    skMat = lin <> predCov <> tr lin + unSym (measCov input)
    kkMat = predCov <> tr lin <> unsafeInv skMat

runKFUpdate
  :: (KnownNat m, KnownNat n)
  => (a -> R n -> L m n) -- ^ Linear measurement operator at a point \(\boldsymbol{h}_i\)
  -> (a -> Sym m)        -- ^ Covariance matrix encoding measurement noise \(R_i\)
  -> a                   -- ^ Dynamical input \(\boldsymbol{u}_{i-1}\)
  -> (R n, Sym n)        -- ^ Current prediction \(({\boldsymbol{x}}_i^\flat, {\boldsymbol{\Sigma}}_i^\flat)\)
  -> R m                 -- ^ New measurement \(\boldsymbol{y}_i\)
  -> (R n, Sym n)        -- ^ Updated prediction \((\hat{\boldsymbol{x}}_i, \hat{\boldsymbol{\Sigma}}_i)\)
runKFUpdate linMeas =
  runEKFUpdate (\inp sys -> linMeas inp sys #> sys) linMeas

-- | Here we combine the prediction and update setps applied to a new
-- measurement, thereby creating a single step of the (extended) Kalman
-- filter.
runEKF
  :: (KnownNat m, KnownNat n)
  => (a -> R n -> R m)   -- ^ System measurement function \(\boldsymbol{h}_i\)
  -> (a -> R n -> L m n) -- ^ Linearization of the measurement at a point
                         -- \(\frac{\partial \boldsymbol{h}_i}{\partial \boldsymbol{x}}\big|_{\boldsymbol{x}}\)
  -> (a -> Sym m)        -- ^ Covariance matrix encoding measurement noise \(R_i\)
  -> (a -> R n -> R n)   -- ^ System evolution function
                         -- \(\boldsymbol{a}_i(\boldsymbol{u}, \boldsymbol{x})\),
                         -- note the order of the control and the
                         -- input
  -> (a -> R n -> Sq n)  -- ^ Linearization of the system evolution at a point  \(\frac{\partial \boldsymbol{a}_i}{\partial \boldsymbol{x}}\big|_{{\boldsymbol{u}}, \, {\boldsymbol{x}}}\)
  -> (a -> Sym n)        -- ^ Covariance matrix encoding system evolution noise \(Q_i\)
  -> a                   -- ^ Dynamical input \(\boldsymbol{u}_{i-1}\)
  -> (R n, Sym n)        -- ^ Current estimate
                         -- \((\hat{\boldsymbol{x}}_{i-1}, \hat{\boldsymbol{\Sigma}}_{i-1})\)
  -> R m                 -- ^ New measurement \(\boldsymbol{y}_i\)
  -> (R n, Sym n)        -- ^ New (filtered) estimate \((\hat{\boldsymbol{x}}_i, \hat{\boldsymbol{\Sigma}}_i)\)
runEKF measure linMeas measCov
  evolve linEvol sysCov
  input estSys newMeas = updatedEstimate
  where
    predictedSystem = runEKFPrediction evolve linEvol sysCov input estSys
    updatedEstimate = runEKFUpdate measure linMeas measCov input predictedSystem newMeas

-- | The ordinary Kalman filter is a special case of the extended Kalman
-- filter, when the state update and measurement transformations are
-- already linear.
runKF
  :: (KnownNat m, KnownNat n)
  => (a -> R n -> L m n) -- ^ Linear measurement operator at a point \(\boldsymbol{h}_i\)
  -> (a -> Sym m)        -- ^ Covariance matrix encoding measurement noise \(R_i\)
  -> (a -> R n -> Sq n)  -- ^ Linear system evolution at a point
                         -- \(\boldsymbol{a}_i(\boldsymbol{u}_{i-1}, \boldsymbol{x}_{i-1})\),
                         -- note the order of the control and the
                         -- input
  -> (a -> Sym n)        -- ^ Covariance matrix encoding system evolution noise \(Q_i\)
  -> a                   -- ^ Dynamical input \(\boldsymbol{u}_{i-1}\)
  -> (R n, Sym n)        -- ^ Current estimate
                         -- \((\hat{\boldsymbol{x}}_{i-1}, \hat{\boldsymbol{\Sigma}}_{i-1})\)
  -> R m                 -- ^ New measurement \(\boldsymbol{y}_i\)
  -> (R n, Sym n)        -- ^ New (filtered) estimate \((\hat{\boldsymbol{x}}_i, \hat{\boldsymbol{\Sigma}}_i)\)
runKF linMeas measCov
  linEvol sysCov
  input estSys newMeas = updatedEstimate
  where
    predictedSystem = runKFPrediction linEvol sysCov input estSys
    updatedEstimate = runKFUpdate linMeas measCov input predictedSystem newMeas

runUKFPrediction
  :: KnownNat n
  => (a -> R n -> R n) -- ^ System evolution function at a point
                       -- \(\boldsymbol{a}_i(\boldsymbol{u}, \boldsymbol{x})\),
                       -- note the order of the control and the
                       -- input
  -> (a -> Sym n)      -- ^ Covariance matrix encoding system evolution noise \(Q_i\)
  -> a                 -- ^ Dynamical input \(\boldsymbol{u}_{i-1}\)
  -> (R n, Sym n)      -- ^ Current estimate \((\hat{\boldsymbol{x}}_{n-1}, \hat{\boldsymbol{\Sigma}}_{n-1})\)
  -> (R n, Sym n)      -- ^ Prediction \(({\boldsymbol{x}_i^\flat}, {\boldsymbol{\Sigma}_i^\flat})\)

runUKFPrediction evolve sysCov input (estMu, estCov) =
  (predMu, predCov)
  where
    predMu  = weightM0 * estMu' +
              sum (map (weightCM *) sigmaPoints')

    predCov = sym $ col (weightC0 * (estMu' - predMu)) <> row (estMu' - predMu) +
              sum (map (\sig ->
                          col (weightCM * (sig - predMu)) <> row (sig - predMu)
                       )
                   sigmaPoints') +
              unSym (sysCov input)

    estMu' = evolve input estMu

    sqCov  = chol estCov
    sqRows = map (* sqrt 3) $ toRows sqCov
    sigmaPoints = map (estMu +) sqRows ++ map (estMu -) sqRows -- 2n points
    sigmaPoints' = map (evolve input) sigmaPoints

    -- hand-tuned weights, more explanation required
    n = fromIntegral $ size estMu
    weightM0 = 1 - n / 3
    weightCM = 1 / 6
    weightC0 = 4 - n / 3 - 3 / n

runUKFUpdate
  :: (KnownNat n, KnownNat m)
  => (a -> R n -> R m) -- ^ Measurement transformation \(\boldsymbol{h}_i\)
  -> (a -> Sym m)      -- ^ Covariance matrix encoding measurement noise \(R_i\)
  -> a                 -- ^ Dynamical input \(\boldsymbol{u}_{i-1}\)
  -> (R n, Sym n)      -- ^ Current prediction \(({\boldsymbol{x}}_i^\flat, {\boldsymbol{\Sigma}}_i^\flat)\)
  -> R m               -- ^ New measurement \(\boldsymbol{y}_i\)
  -> (R n, Sym n)      -- ^ Updated prediction \((\hat{\boldsymbol{x}}_i, \hat{\boldsymbol{\Sigma}}_i)\)

runUKFUpdate measure measCov input (predMu, predCov) newMeas =
  (newMu, newCov)
  where
    newMu  = predMu + kkMat #> (newMeas - upMu)
    newCov = sym $ (unSym predCov) - kkMat <> skMat <> tr kkMat

    predMu' = measure input predMu

    kkMat = ckMat <> unsafeInv skMat
    upMu  = weightM0 * predMu' +
            sum (map (weightCM *) sigmaPoints')

    skMat = (unSym $ measCov input) +
            (col $ weightC0 * (predMu' - upMu)) <> (row $ predMu' - upMu) +
            sum (map (\sig ->
                        (col $ weightCM * (sig - upMu)) <> (row $ sig - upMu)
                     )
                 sigmaPoints')

    ckMat = sum $
            zipWith (\preds meas ->
                       (col $ weightCM' * (preds - predMu)) <> (row $ meas - upMu)
                    )
            sigmaPoints sigmaPoints'

    sqCov   = chol predCov
    sqRows = map (* sqrt 3) $ toRows sqCov
    sigmaPoints = map (predMu +) sqRows ++ map (predMu -) sqRows -- 2n points
    sigmaPoints' = map (measure input) sigmaPoints

    -- hand tuned weights, more exlanation required
    n = fromIntegral $ size predMu
    weightM0 = 1 - n / 3
    weightCM = 1 / 6
    weightCM' = 1 / 6
    weightC0 = 4 - n / 3 - 3 / n

runUKF
  :: (KnownNat m, KnownNat n)
  => (a -> R n -> R m) -- ^ System measurement function \(\boldsymbol{h}_i\)
  -> (a -> Sym m)      -- ^ Covariance matrix encoding measurement noise \(R_i\)
  -> (a -> R n -> R n) -- ^ System evolution function
                       -- \(\boldsymbol{a}_i(\boldsymbol{u}, \boldsymbol{x})\),
                       -- note the order of the control and the
                       -- input
  -> (a -> Sym n)      -- ^ Covariance matrix encoding system evolution noise \(Q_i\)
  -> a                 -- ^ Dynamical input \(\boldsymbol{u}_{i-1}\)
  -> (R n, Sym n)      -- ^ Current estimate
                       -- \((\hat{\boldsymbol{x}}_{i-1}, \hat{\boldsymbol{\Sigma}}_{i-1})\)
  -> R m               -- ^ New measurement \(\boldsymbol{y}_i\)
  -> (R n, Sym n)      -- ^ New (filtered) estimate \((\hat{\boldsymbol{x}}_i, \hat{\boldsymbol{\Sigma}}_i)\)
runUKF measure measCov
  evolve sysCov
  input estSys newMeas = updatedEstimate
  where
    predictedSystem = runUKFPrediction evolve sysCov input estSys
    updatedEstimate = runUKFUpdate measure measCov input predictedSystem newMeas

-- | The Kalman smoothing process (sometimes also called the
-- Rauch-Tung-Striebel smoother or RTSS) is a recursive procedure for
-- improving previous estimates of the system state in light of more
-- recently obtained measurements.  That is, for a given state estimate
-- that was produced (by a Kalman filter) using only the available
-- historical data, we improve the estimate (smoothing it) using
-- information that only became available later, via additional
-- measurements after the estimate was originally made.  The result is an
-- estimate of the state of the system at that moment in its evolution,
-- taking full advantage of hindsight obtained from observations of the
-- system in the future of the current step.
--
-- Consequently, the recursive smoothing procedure progresses backwards
-- through the state evolution, beginning with the most recent
-- observation and updating past observations given the newer
-- information.  The update is made by measuring the deviation between
-- what the system actually did (via observation) and what our best
-- estimate of the system at that time predicted it would do, then
-- adjusting the current estimate in view of this deviation.
runEKS
  :: (KnownNat n)
  => (a -> R n -> R n)  -- ^ System evolution function
                        -- \(\boldsymbol{a}_i(\boldsymbol{u}, \boldsymbol{x})\),
                        -- note the order of the control and the
                        -- input
  -> (a -> R n -> Sq n) -- ^ Linearization of the system evolution at a point  \(\frac{\partial \boldsymbol{a}_i}{\partial \boldsymbol{x}}\big|_{{\boldsymbol{u}}, \, {\boldsymbol{x}}}\)
  -> (a -> Sym n)       -- ^ Covariance matrix encoding system evolution noise \(Q_i\)
  -> a                  -- ^ Dynamical input \(\boldsymbol{u}_{i}\)
  -> (R n, Sym n)       -- ^ Future smoothed estimate
                        -- \((\bar{\boldsymbol{x}}_{i+1}, \bar{\boldsymbol{\Sigma}}_{i+1})\)
  -> (R n, Sym n)       -- ^ Present filtered estimate
                        -- \((\hat{\boldsymbol{x}}_{i}, \hat{\boldsymbol{\Sigma}}_{i})\)
  -> (R n, Sym n)       -- ^ Present smoothed estimate
                        -- \((\bar{\boldsymbol{x}}_{i}, \bar{\boldsymbol{\Sigma}}_{i})\)

runEKS sysEvol linEvol sysCov input (futMu, futCov') (curMu, curCov') =
  (smMu, smCov)
  where
    smMu  = curMu + gkMat #> (futMu - predMu)
    smCov = sym $ curCov + gkMat <> (futCov - predCov) <> tr gkMat

    futCov = unSym futCov'

    curCov = unSym curCov'

    (predMu, predCov') = runEKFPrediction sysEvol linEvol sysCov input (curMu, curCov')
    predCov  = unSym predCov'

    gkMat = curCov <> tr lin <> unsafeInv predCov
    lin   = linEvol input curMu

runKS
  :: (KnownNat n)
  => (a -> R n -> Sq n) -- ^ Linear system evolution at a point
                        -- \(\boldsymbol{a}_i(\boldsymbol{u}, \boldsymbol{x})\),
                        -- note the order of the control and the
                        -- input
  -> (a -> Sym n)       -- ^ Covariance matrix encoding system evolution noise \(Q_i\)
  -> a                  -- ^ Dynamical input \(\boldsymbol{u}_{i}\)
  -> (R n, Sym n)       -- ^ Future smoothed estimate
                        -- \((\bar{\boldsymbol{x}}_{i+1}, \bar{\boldsymbol{\Sigma}}_{i+1})\)
  -> (R n, Sym n)       -- ^ Present filtered estimate
                        -- \((\hat{\boldsymbol{x}}_{i}, \hat{\boldsymbol{\Sigma}}_{i})\)
  -> (R n, Sym n)       -- ^ Present smoothed estimate
                        -- \((\bar{\boldsymbol{x}}_{i}, \bar{\boldsymbol{\Sigma}}_{i})\)
runKS linEvol =
  runEKS (\inp sys -> linEvol inp sys #> sys) linEvol

-- | Unscented Kalman smoother
runUKS
  :: (KnownNat n)
  => (a -> R n -> R n) -- ^ System evolution function
                       -- \(\boldsymbol{a}_i(\boldsymbol{u}, \boldsymbol{x})\),
  -> (a -> Sym n)      -- ^ Covariance matrix encoding system evolution noise \(Q_i\)
  -> a                 -- ^ Dynamical input \(\boldsymbol{u}_{i}\)
  -> (R n, Sym n)      -- ^ Future smoothed estimate
                        -- \((\bar{\boldsymbol{x}}_{i+1}, \bar{\boldsymbol{\Sigma}}_{i+1})\)
  -> (R n, Sym n)      -- ^ Present filtered estimate
                        -- \((\hat{\boldsymbol{x}}_{i}, \hat{\boldsymbol{\Sigma}}_{i})\)
  -> (R n, Sym n)      -- ^ Present smoothed estimate
                       -- \((\bar{\boldsymbol{x}}_{i}, \bar{\boldsymbol{\Sigma}}_{i})\)
runUKS evolve sysCov input (futMu, futCov') (curMu, curCov') =
  (smMu, smCov)
  where
    smMu  = curMu + gkMat #> (futMu - predMu)
    smCov = sym $ curCov + gkMat <> (futCov - predCov) <> tr gkMat

    futCov = unSym futCov'

    curCov = unSym curCov'

    (predMu, predCov') = runUKFPrediction evolve sysCov input (curMu, curCov')
    predCov  = unSym predCov'

    gkMat = dkMat <> unsafeInv predCov

    dkMat = sum $
            zipWith (\pres fut ->
                       (col $ weightCM * (pres - curMu)) <> (row $ fut - predMu)
                    )
            sigmaPoints sigmaPoints'

    sqCov   = chol predCov'
    sqRows = map (* sqrt 3) $ toRows sqCov
    sigmaPoints = map (curMu +) sqRows ++ map (curMu -) sqRows -- 2n points
    sigmaPoints' = map (evolve input) sigmaPoints

    -- hand tuned weights, more exlanation required
    weightCM = 1 / 6

unsafeInv :: KnownNat n => Sq n -> Sq n
unsafeInv m = fromJust $ linSolve m eye

max_stars_repo_name: idontgetoutmuch/Kalman

max_issues_repo_name: f-o-a-m/Kalman

max_forks_repo_name: f-o-a-m/Kalman


text:
{-# LANGUAGE GADTs #-}
module Format.FastTime
  ( FastTime()
  , fastTime
  , ignoreHigh
  ) where

import Data.Time.Clock
import Statistics.Distribution
import Statistics.Distribution.Normal
import Text.Printf

newtype FastTime where
  FastTime :: Double -> FastTime

instance Show FastTime where
  showsPrec _ = showFastTime

fastTime :: NominalDiffTime -> FastTime
fastTime = FastTime . realToFrac

showFastTime :: FastTime -> ShowS
showFastTime (FastTime p) = showString go
  where go | p <= (-1000) = replicate 4 '\8593'
           | p >=      0  = replicate 4 '-'
           | p <=   (-10) = printf "-%3.0f" $ abs p
           | otherwise    = printf "%4.1f" p

ignoreHigh :: FastTime -> String
ignoreHigh a = go $ calc a
  where go p | p <= (-10) = "   "
             | otherwise  = show a

calc :: FastTime -> Double
calc (FastTime t) = undefined

ext: hs

max_stars_repo_name: argiopetech/timer

max_issues_repo_name: argiopetech/timer

max_forks_repo_name: argiopetech/timer


text:
{-# LANGUAGE BangPatterns #-}

module School.FileIO.BinConversion
( binConversion
, decodeList
) where

import Data.ByteString (ByteString, drop, length, take, unpack)
import Numeric.LinearAlgebra (I, R)
import Prelude hiding (drop, length, take)
import School.Types.Decoding (binToInt)
import School.Types.Encoding (Put, putInt, putDouble, runPut)
import School.Types.DataType (DataType(..))
import School.Types.Error (Error)

type Converter = ByteString
              -> Either Error ByteString

floatMsg :: Error
floatMsg = "Illegal conversion from floating point to integral"

precMsg :: Error
precMsg = "Illegal conversion from higher to lower "
       ++ "precision variables"

convertWords :: (Num a)
             => (Int -> a)
             -> (a -> Put)
             -> Converter
convertWords fromInt put = Right
                         . runPut
                         . mapM_ (put . fromInt . fromEnum)
                         . unpack

decodeList :: Int
           -> (a -> b)
           -> (ByteString -> Either Error a)
           -> (ByteString -> Either Error [b])
decodeList s convert decode bytes = reverse <$>
  go s bytes [] where
    go n b acc = if length b < n
                   then return acc
                   else do
                     e <- decode . take n $ b
                     let e' = convert e
                     go n (drop n b) (e':acc)

binConversion :: DataType
              -> DataType
              -> Converter
binConversion INT08B INT32B =
  convertWords (fromIntegral :: Int -> I) putInt
binConversion INT08B DBL64B =
  convertWords (fromIntegral :: Int -> R) putDouble
binConversion INT32B INT08B = const $ Left precMsg
binConversion INT32B DBL64B = \bytes -> do
  !doubles <- decodeList 4 fromIntegral binToInt bytes
  return . runPut $ mapM_ putDouble doubles
binConversion DBL64B INT08B = const $ Left floatMsg
binConversion DBL64B INT32B = const $ Left floatMsg
binConversion _ _ = Right

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
module Main where

import qualified Numeric.LinearAlgebra as LA
import qualified MachineLearning.Types as T
import qualified MachineLearning as ML
import qualified MachineLearning.Classification.Binary as CB
import qualified MachineLearning.Classification.OneVsAll as OVA
import qualified MachineLearning.TerminalProgress as TP
import qualified MachineLearning.Optimization as Opt
import MachineLearning.MultiSvmClassifier
import MachineLearning.SoftmaxClassifier

featuresMapParameter = 2

-- processFeatures :: T.Matrix -> T.Matrix
processFeatures muSigma = ML.addBiasDimension . ML.featureNormalization muSigma . ML.mapFeatures featuresMapParameter

calcAccuracy :: (Model a) => a -> T.Matrix -> T.Vector -> T.Vector -> Double
calcAccuracy m x y theta = OVA.calcAccuracy y yPredicted
  where yPredicted = hypothesis m x theta

main = do
  putStrLn "\n=== Optical Recognition of Handwritten Digits Data Set ===\n"
  -- Step 1. Data loading.
  -- Step 1.1 Training Data loading.
  (x, y) <- pure ML.splitToXY <*> LA.loadMatrix "digits_classification/optdigits.tra"
  -- Step 1.1 Testing Data loading.
  (xTest, yTest) <- pure ML.splitToXY <*> LA.loadMatrix "digits_classification/optdigits.tes"

  -- Step 2. Outputs and features preprocessing.
  let numLabels = 10
      svm = MultiClass $ MultiSvm 1 numLabels
      softmax = MultiClass $ Softmax numLabels
      muSigma = ML.meanStddev (ML.mapFeatures featuresMapParameter x)
      x1 = processFeatures muSigma x
      xTest1 = processFeatures muSigma xTest
      initialTheta = LA.konst 0.1 (numLabels*(LA.cols x1))
  print $ LA.size x1
  print $ LA.size $ x1 LA.? [1..10]
  -- Step 3. Learning.
  putStrLn "Learning Multi SVM model..."
  (thetaSvm, optPathSvm) <- TP.learnWithProgressBar (Opt.minimize (Opt.BFGS2 0.01 0.7) svm 0.0001 5 (CB.L2 2) x1 y) initialTheta 20

  putStrLn "\nLearning MGD Softmax model"
  --(thetaSm, optPathSm) <- TP.learnWithProgressBar (Opt.minimize (Opt.BFGS2 0.01 0.7) softmax 0.000001 5 1 x1 y) initialTheta 20
  (thetaSm, optPathSm) <- TP.learnWithProgressBar (Opt.minimize (Opt.MinibatchGradientDescent 0 1024 0.01) svm 0.000001 200 (CB.L2 1) x1 y) initialTheta 20

  putStrLn "\nLearning GD Softmax model"
  (thetaSmGD, optPathSmGD) <- TP.learnWithProgressBar (Opt.minimize (Opt.GradientDescent 0.01) softmax 0.000001 50 (CB.L2 1) x1 y) initialTheta 20

  -- Step 4. Prediction and checking accuracy
  let accuracyTrainSvm = calcAccuracy svm x1 y thetaSvm
      accuracyTestSvm = calcAccuracy svm xTest1 yTest thetaSvm

  -- Step 5. Printing results.
  putStrLn $ "\nSVM: Number of iterations to learn: " ++ show (LA.rows optPathSvm)

  putStrLn $ "SVM: Accuracy on train set (%): " ++ show (accuracyTrainSvm*100)
  putStrLn $ "SVM: Accuracy on test set (%): " ++ show (accuracyTestSvm*100)

  let accuracyTrainSm = calcAccuracy softmax x1 y thetaSm
      accuracyTestSm = calcAccuracy softmax xTest1 yTest thetaSm

  -- Step 5. Printing results.
  putStrLn $ "\nSoftmax: Number of iterations to learn: " ++ show (LA.rows optPathSm)

  putStrLn $ "Softmax: Accuracy on train set (%): " ++ show (accuracyTrainSm*100)
  putStrLn $ "Softmax: Accuracy on test set (%): " ++ show (accuracyTestSm*100)

  let accuracyTrainSmGD = calcAccuracy softmax x1 y thetaSmGD
      accuracyTestSmGD = calcAccuracy softmax xTest1 yTest thetaSmGD

  -- Step 5. Printing results.
  putStrLn $ "\nGD Softmax: Number of iterations to learn: " ++ show (LA.rows optPathSmGD)

  putStrLn $ "GD Softmax: Accuracy on train set (%): " ++ show (accuracyTrainSmGD*100)
  putStrLn $ "GD Softmax: Accuracy on test set (%): " ++ show (accuracyTestSmGD*100)

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Main where

import Control.Applicative
import Graphics.Rendering.Chart.Backend.Diagrams (toFile)
import qualified Graphics.Rendering.Chart.Easy as Plt
import Kalman
import Numeric.LinearAlgebra
import Sim

fA1 :: (Floating a) => [a] -> [a]
fA1 [x, y, xh, yh, ax, ay] =
  [ x + xh * dt + c * ax * dt * dt,
    y + yh * dt + c * ay * dt * dt,
    xh + ax * dt,
    yh + ay * dt,
    ax,
    ay
  ]
  where
    dt = 0.1
    c = 0.5

fH1 :: (Floating a) => [a] -> [a]
fH1 [x, y, xh, yh, ax, ay] = [x, y]

cannonBall = toFile Plt.def "CannonBall.svg" $ do
  -- KF Initial conditions
  let x0 = [0, 0, 10, 10, 0, -9.8]
  let q = diagl [0.1, 0.1, 1e-3, 1e-3, 1e-3, 1e-3]
  let r = diagl (replicate 2 0.5)
  let s0 = Estimate (vector x0) (diagl (replicate 6 5))

  -- Fire the cannon ball
  let plant = Plant fA1 id fH1 id
  let sim = take 50 $ simulate plant (vector x0) 0
  let measurements = map snd sim
  let trajectory = map fst sim

  -- Run KF/EKF/UKF
  let sys = System fA1 fH1 q r
  let ekfs = scanl (ekf sys) s0 measurements
  let ukfs = scanl (ukf sys) s0 measurements

  -- Plot the results
  let vector2Chart = map ((\[x, y] -> (x, y)) . toList . subVector 0 2)
  let state2Chart xs = vector2Chart $ map sX xs

  Plt.plot (Plt.line "Plant" [vector2Chart trajectory])
  Plt.plot (Plt.line "UKF" [state2Chart ukfs])
  Plt.plot (Plt.line "EKF" [state2Chart ekfs])
  Plt.plot (Plt.points "measurements" (vector2Chart measurements))

main = cannonBall


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  RL.Util
-- Description :  Misc. Utilities
-- Copyright   :  (c) David Banas, 2018
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  David.Banas@target.com
-- Stability   :  experimental
-- Portability :  ?
--
-- Common utilities used throughout the RL package.
-----------------------------------------------------------------------------

{-# OPTIONS_GHC -Wall #-}

{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE TypeApplications #-}

module RL.Util where

import qualified Prelude as P
import Prelude (Show(..))
import Protolude  hiding (show, for)

import GHC.TypeLits

import qualified Data.Vector.Sized   as VS
import Data.Vector.Sized             (Vector)

import Control.Monad.Writer
import Data.Finite
import Data.Finite.Internal
import Data.List                     ((!!))
import Statistics.Distribution       (density)
import Statistics.Distribution.Gamma (gammaDistr)
import Text.Printf

import ConCat.TArr

-- | Convenient abbreviation of: `natVal (Proxy @n)`.
nat :: forall n. KnownNat n => Integer
nat = natVal (Proxy @n)

-- | Convenient abbreviation of: `fromIntegral (nat @n)`.
int :: forall n. KnownNat n => Int
int = fromIntegral (nat @n)

-- | Apply the matrix representation of a two argument function.
--
-- The first argument is assumed to index the rows of the matrix.
appFm
  :: ( HasFin' a,   HasFin' b )
  => Vector (Card a) (Vector (Card b) c)  -- ^ matrix representation of @f(r,c)@
  -> a                                    -- ^ row
  -> b                                    -- ^ column
  -> c
appFm f x y = f `VS.index` toFin x `VS.index` toFin y

-- | Apply the vector representation of a function.
appFv
  :: ( HasFin' x )
  => Vector (Card x) a  -- ^ vector representation of @f(x)@
  -> x                  -- ^ function argument
  -> a
appFv f x = f `VS.index` toFin x

-- | Convert an action-value matrix to a value vector.
--
-- (i.e. - \(Q(s,a) -> V(s)\))
qToV
  :: ( Ord a
     , KnownNat m, KnownNat n, KnownNat k
     , n ~ (k + 1)
     )
  => Vector m (Vector n a)  -- ^ matrix representation of @Q(s,a)@
  -> Vector m a             -- ^ vector representation of @V(s)@
qToV = VS.map VS.maximum

-- | Convert an action-value matrix to a policy vector.
--
-- (i.e. - \(Q(s,a) -> A(s)\))
qToP
  :: ( Ord a, HasFin' act
     , KnownNat m, KnownNat n, KnownNat k
     , n ~ (k + 1), Card act ~ n
     )
  => Vector m (Vector n a)  -- ^ matrix representation of Q(s,a)
  -> Vector m act           -- ^ vector representation of P(s)
qToP = VS.map (unFin . VS.maxIndex)

{----------------------------------------------------------------------
  Misc.
----------------------------------------------------------------------}

-- | To control the formatting of printed floats in output matrices.
newtype Pfloat = Pfloat { unPfloat :: Float}
  deriving (Eq)

instance Show Pfloat where
  show x = printf "%4.1f" (unPfloat x)

-- | To control the formatting of printed doubles in output matrices.
newtype Pdouble = Pdouble { unPdouble :: Double }
  deriving (Eq, Ord)

instance Show Pdouble where
  show x = printf "%4.1f" (unPdouble x)

poisson :: Finite 5 -> Finite 12 -> Float
poisson (Finite lambda) (Finite n') =
  lambda' ^ n * exp (-lambda') / fromIntegral (fact n)
 where lambda' = fromIntegral lambda
       n       = fromIntegral n'

fact :: Int -> Int
fact 0 = 1
fact n = product [1..n]

poissonVals :: VS.Vector 5 (VS.Vector 12 Float)
poissonVals = VS.generate (VS.generate . poisson)

poisson' :: Finite 5 -> Finite 21 -> Float
poisson' n x@(Finite x') =
  if x > 11  -- The Python code enforces this limit. And we're trying
    then 0   -- for an "apples-to-apples" performance comparison.
    else poissonVals `VS.index` n `VS.index` finite x'

-- | Gamma pdf
--
-- Assuming `scale = 1`, `shape` should be: 1 + mean.
gammaPdf :: Double -> Double -> Double -> Double
gammaPdf shape scale = density (gammaDistr shape scale)

-- | Gamma pmf
--
-- Scale assumed to be `1`, so as to match the calling signature of
-- `poisson`.
gamma :: Finite 5 -> Finite 12 -> Double
gamma (Finite expect') (Finite n') =
  0.1 * sum [gammaPdf' (n + x) | x <- [0.1 * m | m <- [-4..5]]]
    where gammaPdf' = gammaPdf (1 + expect) 1
          expect    = fromIntegral expect'
          n         = fromIntegral n'

gammaVals :: VS.Vector 5 (VS.Vector 12 Double)
gammaVals = VS.generate (VS.generate . gamma)

gamma' :: Finite 5 -> Finite 21 -> Double
gamma' n (Finite x') =
  if x' > 11
    then 0
    else gammaVals `VS.index` n `VS.index` finite x'

-- | Monadically search list for first element less than
-- given threshold under the given function, and return the last element
-- if the threshold was never met.
-- Return 'Nothing' if the input list was empty.
withinOnM :: Monad m
          => Double
          -> (a -> m Double)
          -> [a]
          -> m (Maybe a)
withinOnM _   _ [] = return Nothing
withinOnM eps f xs = do
  n <- withinIxM eps $ map f xs
  case n of
    Nothing -> return Nothing
    Just n' -> return $ Just (xs !! n')

-- | Monadically find index of first list element less than or equal to
-- given threshold, or the index of the last element if the threshold
-- was never met.
--
-- A return value of 'Nothing' indicates an empty list was given.
withinIxM :: Monad m
          => Double
          -> [m Double]
          -> m (Maybe Int)
withinIxM _   [] = return Nothing
withinIxM eps xs = withinIxM' 0 xs
 where withinIxM' n []     = return $ Just (n - 1)
       withinIxM' n (y:ys) = do
         y' <- y
         if y' < eps then return (Just n)
                     else withinIxM' (n+1) ys

-- | Return the maximum value of a set, as well as a count of the number
-- of non-zero elements in the set.
--
-- (See documentation for `chooseAndCount` function.)
maxAndNonZero :: (Foldable t, Num a, Ord a) => t a -> Writer [Int] a
maxAndNonZero = chooseAndCount max (/= 0)

-- | Choose a value from the set using the given comparison function,
-- and provide a count of the number of elements in the set meeting the
-- given criteria.
chooseAndCount :: (Foldable t, Num a)
               => (a -> a -> a)  -- ^ choice function
               -> (a -> Bool)    -- ^ counting predicate
               -> t a            -- ^ foldable set of elements to count/compare
               -> Writer [Int] a
chooseAndCount f p xs = do
  let (val, cnt::Int) =
        foldl' ( \ (v, c) x ->
                   ( f v x
                   , if p x
                       then c + 1
                       else c
                   )
               ) (0,0) xs
  tell [cnt]
  return val

for :: [a] -> (a -> b) -> [b]
for = flip map

vsFor :: Vector n a -> (a -> b) -> Vector n b
vsFor = flip VS.map

-- | Mean value of a collection
mean :: (Foldable f, Fractional a) => f a -> a
mean = uncurry (/) . second fromIntegral . foldl' (\ (!s, !n) x -> (s+x, n+1)) (0,0::Integer)

-- | Find the mean square of a list of lists.
arrMeanSqr :: (Functor f, Foldable f, Functor g, Foldable g, Fractional a) => f (g a) -> a
arrMeanSqr = mean . fmap mean . fmap (fmap sqr)

-- | Take the square of a numerical type.
sqr :: Num a => a -> a
sqr x = x * x

-- | Convert any showable type to a string, avoiding the introduction
-- of extra quotation marks when that type is a string to begin with.
toString :: (Show a, Typeable a) => a -> P.String
toString x = fromMaybe (show x) (cast x)

-- | Convert a Bool to a Double.
boolToDouble :: Bool -> Double
boolToDouble True = 1
boolToDouble _    = 0

-- | Take every nth element from a list.
takeEvery :: Int -> [a] -> [a]
takeEvery _ [] = []
takeEvery n xs = P.head xs : (takeEvery n $ drop n xs)

max_stars_repo_name: capn-freako/haskell-rl

max_issues_repo_name: capn-freako/haskell-rl

max_forks_repo_name: capn-freako/haskell-rl


text:
--
-- CNN : CNN filter generator
--
{-# LANGUAGE BangPatterns #-}

module Main (
  main
) where

import Control.DeepSeq
import Control.Monad
--import Data.List
import Data.Time
import Data.List (foldl')
import Debug.Trace
import Numeric.LinearAlgebra
import System.Environment
import System.IO
import Text.Printf

import CNN.ActLayer
import CNN.Algebra
import CNN.ConvLayer
import CNN.FullConnLayer
import CNN.Image
import CNN.Layer
import CNN.LayerType
import CNN.PoolLayer
import CNN.Propagation

import Pool
import Status

usage :: String
usage = "Usage: train <dir> [<image>]"

-- MAIN

main :: IO ()
main = do
  as <- getArgs
  case length as of
    1 -> train (as !! 0)
    2 -> judge (as !! 0) (as !! 1)
    _ -> error usage

train :: String -> IO ()
train dn = do
  putStrLn "Initializing..."
  st <- loadStatus dn
--  st <- loadStatus ""
  tm0 <- getCurrentTime

  let
    --getTeachers = getImages (poolT st) (batchSz st)
    --getTests    = getImages (poolE st) (testSz st)
    getTeachers = getImagesRandomly (poolT st) (batchSz st)
    getTests    = getImagesRandomly (poolE st) (testSz st)
    putF = putStatus tm0 st getTests
    loopFunc = trainLoop' getTeachers putF st
  putStrLn ("Training the model... (#batch:" ++ (show (batchSz st * nclass st)) ++ ")")
  putF 0 (layers st)
  layers' <- foldM' loopFunc (layers st) [1 .. (repeatCt st)]

  putStrLn "Saving status..."
  saveStatus st layers' (repeatCt st * batchSz st)
  putStrLn "Finished!"

judge :: String -> String -> IO ()
judge dn imf = do
  st <- loadStatus dn
  fh <- openFile imf ReadMode
  lst <- hGetContents fh
  forM_ (lines lst) $ \i -> do
    im <- loadImage i
    putStr (i ++ ",")
    let (y, _) = judgeImage (layers st) im
    forM_ (zip [0..] (toList y)) $ \(c, v) -> do
      let res = show c ++ (printf ":%.2f," (v*100))
      putStr res
    putStrLn ""

--
-- FUNCTIONS
--

{- |
foldM' : monadic version of foldl'

  URL: http://stackoverflow.com/questions/8919026/does-haskell-have-foldlm

-}

foldM' :: (Monad m) => ([a] -> b -> m [a]) -> [a] -> [b] -> m [a]
foldM' _ z [] = return z
foldM' f z (x:xs) = do
  z' <- f z x
  z' `seq` foldM' f z' xs

{- |
trainLoop

  IN : func of getting teachers (including batch size and pool)
       sample of evaluation
       output func
       learning rate
       layers
       epoch numbers

  OUT: updated layers

-}

trainLoop' :: (Int -> IO [Trainer]) -> (Int -> [Layer] -> IO ()) -> Status
           -> [Layer] -> Int -> IO [Layer]
trainLoop' getT putF st ls i = do
  --teachers <- getT ((i-1) * batchSz st)
  teachers <- getT 0
  let ls' = updateLayers (learnR st) teachers ls
  if i `mod` (savePt st) == 0
    then putF i ls'
    else return ()
  return ls'

{- |
updateLayers
-}

updateLayers :: Double -> [Trainer] -> [Layer] -> [Layer]
updateLayers lr ts ls = update lr ls (transpose dls)
  where
    rls = tail $ map reverseLayer $ reverse ls    -- fist element isn't used
    dls = map (trainLayers ls rls) ts             -- dls = diff of layers

{-
putStatus

  IN : start time
       status
       start epoch number
       end epoch number
       step size of status output
       epoch
       layers
       sample of evaluation

-}

putStatus :: UTCTime -> Status -> (Int -> IO [Trainer]) -> Int -> [Layer] -> IO ()
putStatus tm0 st getE i ls = do
  --tests <- getE ((i-1) * testSz st)
  tests <- getE 0
  let
    --(rv, rr) = unzip $ evaluate ls tests
    rr = foldl' (evaluateLayers ls) 0.0 tests
  let
    ite = printf "iter = %5d/%d " i (repeatCt st)
    acc = printf "accuracy = %.10f " (rr / fromIntegral (length tests))
  putStr (ite ++ acc)
  tm <- getCurrentTime
  putStrLn ("time = " ++ show (diffUTCTime tm tm0))
  saveStatus st ls (batchSz st * i)
  --mapM_ putOne rs   -- for debug
  where
    putOne :: ([Double], Double) -> IO ()
    putOne (v, r) = putStrLn ("result:" ++ show v ++ ", ratio:" ++ show r)

max_stars_repo_name: eijian/deeplearning

max_issues_repo_name: eijian/deeplearning

max_forks_repo_name: eijian/deeplearning


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE InstanceSigs           #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE NamedFieldPuns         #-}
{-# LANGUAGE QuantifiedConstraints  #-}
{-# LANGUAGE RecordWildCards        #-}
{-# LANGUAGE ScopedTypeVariables    #-}

module Q.MonteCarlo where
import           Control.Monad.State
import           Data.RVar
import           Q.Stochastic.Discretize
import           Q.Stochastic.Process
import           Control.Monad
import           Q.ContingentClaim
import Data.Random
import Q.Time
import Data.Time
import Statistics.Distribution (cumulative)
import Statistics.Distribution.Normal (standard)
import Q.ContingentClaim.Options
import Q.Types

type Path b = [(Time, b)]

-- |Summary type class aggregates all priced values of paths
class (PathPricer p)  => Summary m p | m->p where
  -- | Updates summary with given priced pathes
  sSummarize      :: m -> [p] -> m

  -- | Defines a metric, i.e. calculate distance between 2 summaries
  sNorm           :: m -> m -> Double

-- | Path generator is a stochastic path generator
class PathGenerator m where
  pgMkNew         :: m->IO m
  pgGenerate      :: Integer -> m -> Path b

-- | Path pricer provides a price for given path
class PathPricer m where
  ppPrice :: m -> Path b -> m

type MonteCarlo s a = StateT [(Time, s)] RVar a

-- | Generate a single trajectory stopping at each provided time.
trajectory :: forall a b d. (StochasticProcess a b, Discretize d b) =>
             d        -- ^ Discretization scheme
           -> a        -- ^ The stochastic process
           -> b        -- ^ \(S(0)\)
           -> [Time]   -- ^ Stopping points \(\{t_i\}_i^n \) where \(t_i > 0\)
           -> [RVar b] -- ^ \(dW\)s. One for each stopping point.
           -> RVar [b] -- ^ \(S(0) \cup \{S(t_i)\}_i^n \)
trajectory disc p s0 times dws = reverse <$> evalStateT (onePath times dws) initState' where
  initState' :: [(Time, b)]
  initState' = [(0, s0)]

  onePath :: [Time] -> [RVar b] -> MonteCarlo b [b]
  onePath [] _ = do
    s <- get
    return $ map snd s
  onePath (t1:tn) (dw1:dws) = do
    s <- get
    let t0 = head s
    b <- lift $ pEvolve p disc t0 t1 dw1
    put $ (t1, b) : s
    onePath tn dws

-- | Generate multiple trajectories. See 'trajectory'
trajectories:: forall a b d. (StochasticProcess a b, Discretize d b) =>
             Int        -- ^Num of trajectories
           -> d          -- ^Discretization scheme
           -> a          -- ^The stochastic process
           -> b          -- ^\(S(0)\)
           -> [Time]     -- ^Stopping points \(\{t_i\}_i^n \) where \(t_i > 0\)
           -> [RVar b]   -- ^\(dW\)s. One for each stopping point.
           -> RVar [[b]] -- ^\(S(0) \cup \{S(t_i)\}_i^n \)
trajectories n disc p initState times dws = replicateM n $ trajectory disc p initState times dws

observationTimes :: ContingentClaim a -> [Day]
observationTimes = undefined

class Model a b | a -> b where
  discountFactor :: a -> YearFrac -> YearFrac -> RVar Rate
  evolve   :: a -> YearFrac -> StateT (YearFrac, b) RVar Double

max_stars_repo_name: ghais/lowq

max_issues_repo_name: ghais/lowq

max_forks_repo_name: ghais/lowq


text:
module AsciiString.Prelude
(
  module Exports,
  Product2(..),
  Product3(..),
  Product4(..),
  showBits,
)
where

-- base
-------------------------
import Control.Applicative as Exports
import Control.Arrow as Exports
import Control.Category as Exports
import Control.Concurrent as Exports
import Control.Exception as Exports
import Control.Monad as Exports hiding (mapM_, sequence_, forM_, msum, mapM, sequence, forM)
import Control.Monad.IO.Class as Exports
import Control.Monad.Fix as Exports hiding (fix)
import Control.Monad.ST as Exports
import Data.Bits as Exports
import Data.Bool as Exports
import Data.Char as Exports
import Data.Coerce as Exports
import Data.Complex as Exports
import Data.Data as Exports
import Data.Dynamic as Exports
import Data.Either as Exports
import Data.Fixed as Exports
import Data.Foldable as Exports
import Data.Function as Exports hiding (id, (.))
import Data.Functor as Exports
import Data.Functor.Identity as Exports
import Data.Int as Exports
import Data.IORef as Exports
import Data.Ix as Exports
import Data.List as Exports hiding (sortOn, isSubsequenceOf, uncons, concat, foldr, foldl1, maximum, minimum, product, sum, all, and, any, concatMap, elem, foldl, foldr1, notElem, or, find, maximumBy, minimumBy, mapAccumL, mapAccumR, foldl')
import Data.Maybe as Exports
import Data.Monoid as Exports hiding (Last(..), First(..), (<>))
import Data.Ord as Exports
import Data.Proxy as Exports
import Data.Ratio as Exports
import Data.Semigroup as Exports
import Data.STRef as Exports
import Data.String as Exports
import Data.Traversable as Exports
import Data.Tuple as Exports
import Data.Unique as Exports
import Data.Version as Exports
import Data.Word as Exports
import Debug.Trace as Exports
import Foreign as Exports hiding (void)
import GHC.Conc as Exports hiding (withMVar, threadWaitWriteSTM, threadWaitWrite, threadWaitReadSTM, threadWaitRead)
import GHC.Exts as Exports (lazy, inline, sortWith, groupWith)
import GHC.Generics as Exports (Generic)
import GHC.IO.Exception as Exports
import Numeric as Exports
import Prelude as Exports hiding (concat, foldr, mapM_, sequence_, foldl1, maximum, minimum, product, sum, all, and, any, concatMap, elem, foldl, foldr1, notElem, or, mapM, sequence, id, (.))
import System.Environment as Exports
import System.Exit as Exports
import System.IO as Exports
import System.IO.Error as Exports
import System.IO.Unsafe as Exports
import System.Mem as Exports
import System.Mem.StableName as Exports
import System.Timeout as Exports
import Text.ParserCombinators.ReadP as Exports (ReadP, ReadS, readP_to_S, readS_to_P)
import Text.ParserCombinators.ReadPrec as Exports (ReadPrec, readPrec_to_P, readP_to_Prec, readPrec_to_S, readS_to_Prec)
import Text.Printf as Exports (printf, hPrintf)
import Text.Read as Exports (Read(..), readMaybe, readEither)
import Unsafe.Coerce as Exports

-- bytestring
-------------------------
import Data.ByteString as Exports (ByteString)
import Data.ByteString.Short as Exports (ShortByteString)

-- primitive
-------------------------
import Data.Primitive.PrimArray as Exports hiding (sizeOf, alignment)

-- cereal
-------------------------
import Data.Serialize as Exports

-- foldl
-------------------------
import Control.Foldl as Exports (Fold(..), FoldM(..))

-- deferred-folds
-------------------------
import DeferredFolds.Unfoldl as Exports (Unfoldl(..))
import DeferredFolds.UnfoldlM as Exports (UnfoldlM(..))
import DeferredFolds.Unfoldr as Exports (Unfoldr(..))

-- hashable
-------------------------
import Data.Hashable as Exports

-- deepseq
-------------------------
import Control.DeepSeq as Exports

data Product2 a b = Product2 !a !b

data Product3 a b c = Product3 !a !b !c

data Product4 a b c d = Product4 !a !b !c !d

showBits :: forall a. (Integral a, Show a, FiniteBits a) => a -> String
showBits = padL (finiteBitSize (undefined :: a)) '0' . ($ "") . showIntAtBase 2 intToDigit

padL :: Int -> Char -> String -> String
padL n c s
    | length s < n  = s ++ replicate (n - length s) c
    | otherwise     = s

max_stars_repo_name: metrix-ai/ascii-string

max_issues_repo_name: metrix-ai/ascii-string

max_forks_repo_name: metrix-ai/ascii-string


text:
{-# LANGUAGE FlexibleContexts           #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE RecordWildCards            #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE StandaloneDeriving         #-}

module Main (main)
where

import           Phase
import           Random

import           Control.Monad                     as CM
import qualified Control.Monad.Log                 as ML
import           Control.Monad.Random
import qualified Control.Monad.State.Strict        as MS
import           Data.Complex
import qualified Data.List                         as L
import qualified Data.Massiv.Array                 as A
import qualified Data.Massiv.Array.Manifest.Vector as AMV
import qualified Data.MultiSet                     as M
import qualified Data.Vector.Unboxed               as V
import           Numeric
import           Prelude                           as P
import           System.Console.ANSI
import           System.Environment
import           System.TimeIt

-- | a vector of multiplication phases; should be kept sorted
type PhaseVector = A.Array A.U A.Ix1 Phase

-- | the massiv computation strategy to use throughout
strat :: A.Comp
strat = A.Par

fromVector :: Vector Phase -> PhaseVector
fromVector vec = AMV.fromVector' strat (A.Sz1 (V.length vec)) vec

-- | a tuple of phase vectors to be collimated
type PVTuple = (PhaseVector, PhaseVector)

-- what LoggingT couldn't
deriving instance (MonadRandom m) => MonadRandom (ML.LoggingT message m)

-- | the base-2 logarithm
log2 :: Integral a => a -> Double
log2 = logBase (2 :: Double) . fromIntegral

-- | for convenience
log2rnd :: Integral a => a -> Int
log2rnd  = round . log2

dot :: Num a => [a] -> [a] -> a
dot [] _          = 0
dot _ []          = 0
dot (x:xs) (y:ys) = x*y + dot xs ys

---------- COLLIMATION ----------

data CollimateInfo = SI { s     :: !Phase, lReq :: !Int,
                          lReal :: !Int,   kept :: !Bool }

instance Show CollimateInfo where
  show SI{..} =
    let lRatio = fromIntegral lReal / fromIntegral lReq :: Double
    in "(~2^" P.++ show (log2rnd s) P.++
       "," P.++ show lReq P.++
       "," P.++ show lReal P.++
       "," P.++ showFFloat (Just 2) lRatio
       (if not kept then " DISCARDED" else "") P.++ ")"

collimateInfoHandler :: MonadIO m
                     => Int     -- | number of bottom sieve layers to suppress
                     -> [Phase] -- | list of intervals sizes for sieve
                     -> CollimateInfo
                     -> m ()
collimateInfoHandler skip ss =
  let len = P.length ss in \si@SI{..} ->
    let i' = L.elemIndex s ss
    in maybe (return ())
       (\i -> when (i < len - skip) $ liftIO $ do
           setSGR [SetColor Foreground (if i < 7 then Vivid else Dull)
                   (toEnum (i `mod` 7 + 1))]
           putStrLn $ P.concat (P.replicate i "  ") P.++ show si
           setSGR [Reset]) i'

-- | all the mod-@n@ subset-sums of a vector. Results are not
-- necessarily sorted!
subsetSums :: Vector Phase -> Vector Phase
subsetSums v =
  V.generate (2^V.length v) (go 0 0) where
  go i acc s
    | i >= V.length v = acc
    | otherwise = go (i+1) (acc + if odd s then v V.! i else 0) (s `div` 2)

-- | choose a random phase multiplier from the given phase vector
randomElt :: MonadRandom m => PhaseVector -> m Phase
randomElt v =  (v A.!) <$> getRandomR (0, A.elemsCount v - 1)

-- | the first index (in a sorted vector) of an element that is /at
-- least/ @a@, or 'length v' if none exists
findAtLeast :: PhaseVector -> Phase -> Int
findAtLeast v a = go 0 (A.elemsCount v) where
  go i j | i == j         = i
         | (v A.! i) >= a = i
         -- don't use i+j `div` 2, to avoid overflow on big arrays
         | otherwise      = let k = i + ((j - i) `div` 2)
                            in if (v A.! k) >= a then go i k else go (k+1) j

-- | collimate phase vectors
collimate :: (MonadRandom m)
          => Phase              -- | desired interval size S (upper bound)
          -> PVTuple            -- | phase vectors to collimate
          -> m PhaseVector
collimate s (v1,v2) = do
  b1 <- randomElt v1
  b2 <- randomElt v2
  let q = (b1 + b2) `div` s
      qs = q*s
      start i1 = (i1, findAtLeast v2 (qs - (v1 A.! i1)))
      gen (i1,i2)
        | i1 >= A.elemsCount v1 = Nothing
        | i2 >= A.elemsCount v2 = gen $ start (i1+1)
        | v < qs + s            = Just (v `mod` s, (i1, i2+1))
        | otherwise             = gen $ start (i1+1)
        where v = (v1 A.! i1) + (v2 A.! i2)
      -- Create a Vector because massiv doesn't supported
      -- unknown-length unfoldr.
  return $ A.quicksort $ fromVector $ V.unfoldr gen (start 0)

---------- SIEVE ----------

data SieveState = SS { numQueries :: !Int, maxLength    :: !Int,
                       numNodes   :: !Int, numDiscarded :: !Int}
  deriving (Show, Eq)

newSieveState :: SieveState
newSieveState = SS 0 0 0 0

addQueries :: Int -> SieveState -> SieveState
addQueries i ss@SS{..} = ss { numQueries = numQueries + i }

updateLength :: Int -> SieveState -> SieveState
updateLength l ss@SS{..} = ss { maxLength = max l maxLength }

incrementNodes :: SieveState -> SieveState
incrementNodes ss@SS{..} = ss { numNodes = succ numNodes }

incrementDiscarded :: SieveState -> SieveState
incrementDiscarded ss@SS{..} = ss { numDiscarded = succ numDiscarded }

-- | Collimation sieve.
sieve :: (MonadRandom m, ML.MonadLog CollimateInfo m,
          MS.MonadState SieveState m)
      => Phase   -- | group order N
      -> Double  -- | threshold factor for long-enough phase vectors
      -> [Phase] -- | increasing interval sizes S_i for the i'th level
                 -- of the sieve; the final one should equal N
      -> Int     -- | desired phase vector length L
      -> m PhaseVector
sieve = sieve' True where       -- don't discard final output
  sieve' _ _ _ [] _ = error "sieve: empty list of interval sizes"
  sieve' alwaysKeep n threshold ss@(s:ss') l
    | s >= n = do
        let logl = log2rnd l
        v <- (A.quicksort . fromVector . V.map (`mod` n)) . subsetSums <$>
          -- create a Vector because massiv doesn't support monadic replicate
             V.replicateM logl (getRandomR (0,n-1))
        ML.logMessage (SI s l (A.elemsCount v) True)
        MS.modify'   incrementNodes
        MS.modify' $ addQueries logl
        MS.modify' $ updateLength (A.elemsCount v)
        return v
    | otherwise =
        let s' = P.head ss'
            z  = 1.5 * fromIntegral l * fromIntegral s' / fromIntegral s :: Double
            l' = ceiling $ sqrt z
        in do v1   <- sieve' False n threshold ss' l'
              v2   <- sieve' False n threshold ss' $
                      ceiling (z / fromIntegral (A.elemsCount v1))
              v    <- collimate s (v1,v2)
              let vlen = A.elemsCount v
                  -- check if vector is long enough to be useful
                  keep = alwaysKeep ||
                         fromIntegral vlen / fromIntegral l >= threshold
              ML.logMessage (SI s l vlen keep)
              MS.modify' incrementNodes
              if keep
                then MS.modify' (updateLength vlen) >> return v
                else MS.modify' incrementDiscarded >> sieve' False n threshold ss l

-- | histogram of (n, # multipliers appearing n times) in a phase
-- vector for a (final, small) interval [S]
histogram :: Int -> Vector Int -> [(Int,Int)]
histogram s v =
  let occurs = fmap snd $ M.toOccurList $ M.fromList $ V.toList v
      num    = P.length occurs
      hist   = L.sortOn fst $ M.toOccurList $ M.fromList occurs
  in if s == num then hist else (0, s - num) : hist

chi :: Double -> Complex Double
chi = cis . (2.0 * pi *)

square :: Num a => a -> a
square x = x*x

-- | partition a phase vector into one consisting of its unique phase
-- multipliers, and one containing all the leftovers.
puncture :: PhaseVector -> (PhaseVector, PhaseVector)
puncture v =
  let v' = AMV.toVector v
      u' = V.uniq v'
      gen i
        | i >= V.length v' = Nothing
        | v' V.! (i-1) == v' V.! i = Just (v' V.! i, i+1)
        | otherwise = gen $ i+1
  in (fromVector u', fromVector $ V.unfoldr gen 1)

unfold :: (b -> (a,b)) -> b -> [a]
unfold f b = let (a,b') = f b in a : unfold f b'

-- | probability assigned to \( w \), given by \( \theta = s/N - w/T
-- \), by the \( T \)-dimensional inverse-QFT of the given punctured
-- phase vector on \( [T] \).
probTheta :: Int                -- | range \( [T] \) of phase vector
          -> PhaseVector        -- | uniquified phase vector
          -> Double             -- | \( \theta \)
          -> Double
probTheta t v theta =
  square (magnitude $ A.sum $ A.map (chi . (* theta) . fromIntegral) v) /
  (fromIntegral t * fromIntegral (A.elemsCount v))

-- | probability that the \( T \)-dimensional inverse-QFT on the given
-- punctured phase vector on \( [T] \) will output some \( w \in
-- \lfloor sT/N \rceil - \{ - \lfloor z/2 \rfloor, \ldots, \lfloor
-- (z-1)/2 \rfloor \} \), given the (essentially uniform) shift \(
-- sT/N \bmod 1 \in [-1/2,1/2) \).  Assumes that (\ z \) is positive
-- and small enough not to double-count probabilities.
probClosest :: Int               -- | \( T \)
            -> Double            -- | \( sT/N \bmod 1 \in [-1/2,1/2) \)
            -> PhaseVector       -- | uniquified phase vector
            -> Int               -- | number \( z \) of closest to count
            -> Double
probClosest t shift v z =
  let thetas = A.makeArrayR A.D strat (A.Sz1 z)
               (\(A.Ix1 i) -> (shift + fromIntegral (i - z `div` 2))
                              / fromIntegral t)
  in A.sum $ A.map (probTheta t v) thetas

main :: IO ()
main = do
  args <- getArgs
  let (params,rest) = splitAt 3 args
      [logn :: Int, logl, logs] = read <$> params
      threshold = if null rest then 0.25 :: Double else read (head rest)
      n     = if logn > 0
              then 2^logn
                   -- from https://eprint.iacr.org/2019/498.pdf
              else 3 * 37 * 1407181 * 51593604295295867744293584889 *
                   31599414504681995853008278745587832204909
      l     = 2^logl
      s     = 2^logs
      ss    = P.takeWhile (< n)
              (iterate (\t -> (2 * t * fromIntegral l) `div` 3) $ fromIntegral s)
              P.++ [n]

  -- print schedule of log interval sizes
  putStrLn $ "log S's = " P.++ show (log2rnd <$> ss)

  -- run the sieve
  (time, (v, sieveState@SS{..})) <- timeItT $
    evalCryptoRandIO
    (flip ML.runLoggingT (collimateInfoHandler (min 5 (P.length ss - 1)) ss)
     (flip MS.runStateT newSieveState (sieve n threshold ss l)))

  -- print results

  -- parameters again, for convenience
  putStrLn $ "\n[log N, log L, log S_0] = " P.++
    show [logn, logl, logs]

  putStrLn $ "threshold = " P.++ showFFloat (Just 2) threshold ""

  -- log interval sizes again, for convenience
  putStrLn $ "\nlog S's = " P.++ show (log2rnd <$> ss)

  putStrLn $ "\nSieve summary = " P.++ show sieveState

  -- compute number of queries according to model
  let delta = fromIntegral numDiscarded / fromIntegral numNodes :: Double
      depth = P.length ss - 1
      l' = sqrt $ 1.5 * fromIntegral l * fromIntegral n /
           fromIntegral (ss P.!! (P.length ss - 2)) :: Double
      modelQueries = (2.0 / (1-delta))^depth * logBase 2 l'

  putStrLn $ "\nProbability of discarding = " P.++ showFFloat (Just 4) delta ""

  putStrLn $ (P.++) "\nNumber of queries actual/modeled = " $ (P.++)
    (show numQueries)                $ (P.++) "/" $
    showFFloat (Just 1) modelQueries $ (P.++) " ~= " $
    showFFloat (Just 2) (fromIntegral numQueries / modelQueries) ""

  -- compute probability of obtaining a regular state
  let hist = histogram s $ V.map fromIntegral $ AMV.toVector v
      num  = s * fst (P.head hist)
      den  = A.elemsCount v
      probRegular = fromIntegral num / fromIntegral den :: Double

  putStrLn $ "\nProbability of obtaining a regular state = " P.++
    show num P.++ "/" P.++ show den P.++ " ~= " P.++
    showFFloat (Just 3) probRegular ""

  -- compute probabilities of getting punctured regular states, and close
  shift :: Double <- getRandomR (-0.5, 0.5)
  let punctureds = takeWhile ((> 2^(logs - 6)) . A.elemsCount) $
                   unfold puncture v
      puncturedProbs :: [Double] =
        (/ fromIntegral den) . fromIntegral . A.elemsCount <$> punctureds
      closest = [1,2,4,8]       -- how many of the closest w to check
      closestProbs = map (<$> closest) (probClosest s shift <$> punctureds)
      totalClosestProbs = dot puncturedProbs <$> L.transpose closestProbs

  putStrLn $ (P.++) "\nProbability of obtaining each punctured state (1st, 2nd, ... attempt) =\n" $
    L.intercalate ", " $ mapM (showFFloat (Just 3)) puncturedProbs ""

  putStrLn $ (P.++) "\nProbability bounds for each punctured state for " $
    (P.++) (show closest) $ (P.++) " closest w = \n" $
    L.intercalate "\n" $ L.intercalate ", " <$>
    (mapM . mapM) (showFFloat (Just 3)) closestProbs ""

  putStrLn $ (P.++) "\nTotal probability bounds for " $
    (P.++) (show closest) $ (P.++) " closest w = \n" $
    L.intercalate ", " $ mapM (showFFloat (Just 3)) totalClosestProbs ""

  putStrLn $ (P.++) "\nLaTeX table row:\n" $
    showFFloat (Just 1) (log2 n)                          $ (P.++) " & " $
    showFFloat (Just 1) (log2 numQueries)                 $ (P.++) " & " $
    showFFloat (Just 1) (logBase 2 modelQueries)          $ (P.++) " & " $
    showFFloat (Just 1) (log2 maxLength)                  $ (P.++) " & " $
    (P.++) (show logl)                                    $ (P.++) " & " $
    (P.++) (show logs)                                    $ (P.++) " & " $
    showFFloat (Just 0) (probRegular * 100)               $ (P.++) " & " $
    showFFloat (Just 1) (fromIntegral logs * probRegular) $ (P.++) " & " $
    showFFloat (Just 2) threshold                         $ (P.++) " & " $
    showFFloat (Just 1) (delta * 100)                     $ (P.++) " & " $
    (P.++) (show depth)                                   $ (P.++) " & " $
    showFFloat (Just 1) (time / 3600) " \\\\"

max_stars_repo_name: cpeikert/CollimationSieve

max_issues_repo_name: cpeikert/CollimationSieve

max_forks_repo_name: cpeikert/CollimationSieve


text:
{-# LANGUAGE ScopedTypeVariables #-}
module Main (main) where

import Data.Complex (Complex)
import Data.Int (Int8, Int16, Int32, Int64)
import Data.List.NonEmpty (NonEmpty)
import Data.Proxy (Proxy (..))
import Data.Typeable (Typeable, typeRep)
import Data.Word (Word8, Word16, Word32, Word64)
import Numeric.Natural (Natural)
import Test.QuickCheck (Arbitrary (..), Property, counterexample, label, (===),
                        sized, chooseInt, vectorOf)
import Test.QuickCheck.Instances ()
import Test.Tasty (defaultMain, testGroup, TestTree)
import Test.Tasty.QuickCheck (testProperty)

import qualified Data.List as L

import Data.Discrimination
import Utils

main :: IO ()
main = defaultMain $ testGroup "discrimination"
  [ testGroup "examples"
    [ testGroup "nub"
      [ testProperty "List.nub" $
        let prop :: [Word64] -> Property
            prop xs = L.nub xs === nub xs
        in prop

      , testProperty "ordNub" $
        let prop :: [Word64] -> Property
            prop xs = ordNub xs === nub xs
        in prop

      , testProperty "hashNub" $
        let prop :: [Word64] -> Property
            prop xs = hashNub xs === nub xs
        in prop
      ]

    , testGroup "sort"
      [ testProperty "List.sort" $
        let prop :: [Word64] -> Property
            prop xs = L.sort xs === sort xs
        in prop

      , testProperty "introsort" $
        -- for Word64 unstable sort works too
        let prop :: [Word64] -> Property
            prop xs = introsort xs === sort xs
        in prop

      , testProperty "mergesort" $
        let prop :: [Word64] -> Property
            prop xs = mergesort xs === sort xs
        in prop
      ]
    ]

  , testGroup "Grouping"
    [ testGrouping (Proxy :: Proxy ())
    , testGrouping (Proxy :: Proxy Int)
    , testGrouping (Proxy :: Proxy Int8)
    , testGrouping (Proxy :: Proxy Int16)
    , testGrouping (Proxy :: Proxy Int32)
    , testGrouping (Proxy :: Proxy Int64)
    , testGrouping (Proxy :: Proxy Word)
    , testGrouping (Proxy :: Proxy Word8)
    , testGrouping (Proxy :: Proxy Word16)
    , testGrouping (Proxy :: Proxy Word32)
    , testGrouping (Proxy :: Proxy Word64)
    , testGrouping (Proxy :: Proxy Bool)
    , testGrouping (Proxy :: Proxy Ordering)
    , testGrouping (Proxy :: Proxy (Word8,Word8))
    , testGrouping (Proxy :: Proxy (Word8,Word8,Word8))
    , testGrouping (Proxy :: Proxy (Word8,Word8,Word8,Word8))
    , testGrouping (Proxy :: Proxy Rational)
    , testGrouping (Proxy :: Proxy (Complex Word8))
    , testGrouping (Proxy :: Proxy (Maybe Word8))
    , testGrouping (Proxy :: Proxy (Either Word8 Word8))
    , testGrouping (Proxy :: Proxy Char)
    , testGrouping (Proxy :: Proxy String)
    , testGrouping (Proxy :: Proxy (NonEmpty Int))
    , testGrouping (Proxy :: Proxy Natural)
    , testGrouping (Proxy :: Proxy Integer)

    , testGrouping' listToNatural
    , testGrouping' listToInteger
    ]

  , testGroup "Sorting"
    [ testSorting (Proxy :: Proxy ())
    , testSorting (Proxy :: Proxy Int)
    , testSorting (Proxy :: Proxy Int8)
    , testSorting (Proxy :: Proxy Int16)
    , testSorting (Proxy :: Proxy Int32)
    , testSorting (Proxy :: Proxy Int64)
    , testSorting (Proxy :: Proxy Word)
    , testSorting (Proxy :: Proxy Word8)
    , testSorting (Proxy :: Proxy Word16)
    , testSorting (Proxy :: Proxy Word32)
    , testSorting (Proxy :: Proxy Word64)
    , testSorting (Proxy :: Proxy Bool)
    , testSorting (Proxy :: Proxy Ordering)
    , testSorting (Proxy :: Proxy (Word8,Word8))
    , testSorting (Proxy :: Proxy (Word8,Word8,Word8))
    , testSorting (Proxy :: Proxy (Word8,Word8,Word8,Word8))
    , testSorting (Proxy :: Proxy (Maybe Word8))
    , testSorting (Proxy :: Proxy (Either Word8 Word8))
    , testSorting (Proxy :: Proxy Char)
    , testSorting (Proxy :: Proxy String)
    , testSorting (Proxy :: Proxy (NonEmpty Int))
    , testSorting (Proxy :: Proxy Natural)
    , testSorting (Proxy :: Proxy Integer)

    , testSorting' listToNatural
    , testSorting' listToInteger
    ]
  ]

listToNatural :: SmallList Word64 -> Natural
listToNatural = L.foldl' (\x y -> x * 2 ^ (64 :: Int) + fromIntegral y) 0 . getSmallList

listToInteger :: SmallList Int64 -> Integer
listToInteger = L.foldl' (\x y -> x * 2 ^ (64 :: Int) + fromIntegral y) 0 . getSmallList

newtype SmallList a = SmallList { getSmallList :: [a] } deriving (Eq, Show)

instance Arbitrary a => Arbitrary (SmallList a) where
    arbitrary = sized $ \n -> do
        m <- chooseInt (0, min 10 n)
        SmallList <$> vectorOf m arbitrary

    shrink = fmap SmallList . shrink . getSmallList

testGrouping
  :: forall a. (Grouping a, Typeable a, Arbitrary a, Eq a, Show a)
  => Proxy a
  -> TestTree
testGrouping _ = testGrouping' (id :: a -> a)

testGrouping'
  :: forall a b. (Grouping b, Typeable a, Typeable b, Arbitrary a, Eq b, Show a, Show b)
  => (a -> b)
  -> TestTree
testGrouping' f = testGroup name
    [ testProperty "groupingEq" prop_eq
    , testProperty "nub"        prop_nub
    ]
  where
    tra = typeRep (Proxy :: Proxy a)
    trb = typeRep (Proxy :: Proxy b)
    name = if tra == trb then show tra else show trb ++ " from " ++ show tra

    prop_eq :: a -> a -> Property
    prop_eq x' y' =
        counterexample (show (x,y)) $
        label (show lhs) $
        lhs === groupingEq x y
      where
        x = f x'
        y = f y'
        lhs = x == y

    prop_nub :: [a] -> Property
    prop_nub xs' = L.nub xs === nub xs
      where
        xs = take 100 (map f xs')

testSorting
  :: forall a. (Sorting a, Typeable a, Arbitrary a, Ord a, Show a)
  => Proxy a
  -> TestTree
testSorting _ = testSorting' (id :: a -> a)

testSorting'
  :: forall a b. (Sorting b, Typeable a, Typeable b, Arbitrary a, Ord b, Show a, Show b)
  => (a -> b)
  -> TestTree
testSorting' f = testGroup name
    [ testProperty "sortingCompare" prop_cmp
    , testProperty "sort"           prop_sort
    ]
  where
    tra = typeRep (Proxy :: Proxy a)
    trb = typeRep (Proxy :: Proxy b)
    name = if tra == trb then show tra else show trb ++ " from " ++ show tra

    prop_cmp :: a -> a -> Property
    prop_cmp x' y' =
        counterexample (show (x,y)) $
        label (show lhs) $
        lhs === sortingCompare x y
      where
        x = f x'
        y = f y'
        lhs = compare x y

    prop_sort :: [a] -> Property
    prop_sort xs' = L.sort xs === sort xs
      where
        xs = map f xs'

max_stars_repo_name: bacchanalia/discrimination

max_issues_repo_name: bacchanalia/discrimination

max_forks_repo_name: bacchanalia/discrimination


text:
{-# LANGUAGE FlexibleContexts #-}
module Filter.Pinwheel
  ( module Filter.Utils
  , PinwheelParams(..)
  , pinwheel
  , PinwheelType(..)
  , pinwheelFunc
  , pinwheelFilter
  , convolvePinwheel
  , frequencyDomainMultiply
  ) where

import           Data.Array.Repa     as R
import           Data.Complex
import           Data.List           as L (length, reverse)
import           Data.Vector.Generic as VG (convert)
import           DFT.Plan
import           Filter.Utils
import           Utils.Coordinates

data PinwheelParams = PinwheelParams
  { pinwheelParamsRows         :: !Int
  , pinwheelParamsCols         :: !Int
  , pinwheelParamsAlpha        :: !Double
  , pinwheelParamsRadialPeriod :: !Double
  , pinwheelParamsAngularFreqs :: ![Double]
  , pinwheelParamsRadialFreqs  :: ![Double]
  }

data PinwheelType
  = Pinwheel
  | PinwheelHollow0 !Double
  | PinwheelHollow1 !Double
  deriving (Read, Show)

{-# INLINE rFunc #-}
rFunc :: Int -> Int -> Double
rFunc x y = sqrt . fromIntegral $ x ^ (2 :: Int) + y ^ (2 :: Int)

{-# INLINE thetaFunc #-}
thetaFunc :: Int -> Int -> Double
thetaFunc x y = angleFunctionRad (fromIntegral x) (fromIntegral y)

{-# INLINE pinwheel #-}
pinwheel :: Double -> Double -> Double -> Double -> Int -> Int -> Complex Double
pinwheel af rf maxR alpha x y
  | r == 0 && af == 0 && rf == 0 = 1
  | r == 0 = 0
  | otherwise =
    (r :+ 0) ** (alpha :+ rf * 2 * pi / log maxR) *
    exp (0 :+ af * (thetaFunc x y))
  where
    r = rFunc x y

{-# INCLUDE pinwheelFunc #-}
pinwheelFunc ::
     PinwheelType
  -> Double
  -> Double
  -> Double
  -> Double
  -> Int
  -> Int
  -> Complex Double
pinwheelFunc Pinwheel af rf maxR alpha x y
  | r == 0 && af == 0 && rf == 0 = 1
  | r == 0 = 0
  | otherwise =
    (r :+ 0) ** (alpha :+ rf * 2 * pi / logMaxR) *
    exp (0 :+ af * (thetaFunc x y))
  where
    r = rFunc x y
    logMaxR =
      if maxR == 1
        then 1
        else log maxR
pinwheelFunc (PinwheelHollow0 radius) af rf maxR alpha x y
  | r < radius = 0.0
  | otherwise =
    (r :+ 0) ** (alpha :+ rf * 2 * pi / logMaxR) *
    exp (0 :+ af * (thetaFunc x y))
  where
    r = rFunc x y
    logMaxR =
      if maxR == 1
        then 1
        else log maxR
pinwheelFunc (PinwheelHollow1 radius) af rf maxR alpha x y
  | r == 0 && af == 0 && rf == 0 = 1
  | r < radius = 0.0
  | otherwise =
    (r :+ 0) ** (alpha :+ rf * 2 * pi / logMaxR) *
    exp (0 :+ af * (thetaFunc x y))
  where
    r = rFunc x y
    logMaxR =
      if maxR == 1
        then 1
        else log maxR

{-# INLINE pinwheelPI #-}
pinwheelPI :: Double -> Double -> Double -> Double -> Int -> Int -> Complex Double
pinwheelPI af rf maxR alpha x y
  | r == 0 = 0
  | otherwise =
    (((r) :+ 0) ** (alpha :+ (rf * 2 * pi / (log maxR)))) *
    exp (0 :+ ((af) * theta))
  where
    r = (sqrt . fromIntegral $ x ^ (2 :: Int) + y ^ (2 :: Int))
    theta = pi  + angleFunctionRad (fromIntegral x) (fromIntegral y)

pinwheelFilter ::
     DFTPlan
  -> PinwheelParams
  -> IO (DFTPlan, Array U DIM4 (Complex Double), Array U DIM4 (Complex Double))
pinwheelFilter plan (PinwheelParams rows cols alpha radialPeriod angularFreqs radialFreqs) = do
  let numAngularFreq = L.length angularFreqs
      numRadialFreq = L.length radialFreqs
      filterVec =
        VG.convert . toUnboxed . computeS . makeFilter2D $
        R.traverse2
          (fromListUnboxed (Z :. L.length angularFreqs) angularFreqs)
          (fromListUnboxed (Z :. L.length radialFreqs) radialFreqs)
          (\(Z :. numAngularFreq') (Z :. numRadialFreq') ->
             (Z :. numAngularFreq' :. numRadialFreq' :. cols :. rows))
          (\f1 f2 (Z :. af :. rf :. i :. j) ->
             pinwheel
               (f1 (Z :. af))
               (f2 (Z :. rf))
               radialPeriod
               alpha
               (i - center cols)
               (j - center rows))
      filterPIVec =
        VG.convert . toUnboxed . computeS . makeFilter2D $
        R.traverse2
          (fromListUnboxed (Z :. L.length angularFreqs) angularFreqs)
          (fromListUnboxed (Z :. L.length radialFreqs) radialFreqs)
          (\(Z :. numAngularFreq') (Z :. numRadialFreq') ->
             (Z :. numAngularFreq' :. numRadialFreq' :. cols :. rows))
          (\f1 f2 (Z :. af :. rf :. i :. j) ->
             pinwheelPI
               (f1 (Z :. af))
               (f2 (Z :. rf))
               radialPeriod
               alpha
               (i - div cols 2)
               (j - div rows 2))
      forwardPlanID =
        DFTPlanID DFT1DG [numAngularFreq, numRadialFreq, cols, rows] [2, 3]
      backwardPlanID =
        DFTPlanID IDFT1DG [numAngularFreq, numRadialFreq, cols, rows] [2, 3]
  lock <- getFFTWLock
  (plan1, filterVecF) <-
    dft1dGPlan
      lock
      plan
      [numAngularFreq, numRadialFreq, cols, rows]
      [2, 3]
      filterVec
  (newPlan, _) <-
    idft1dGPlan
      lock
      plan1
      [numAngularFreq, numRadialFreq, cols, rows]
      [2, 3]
      filterVecF
  filterPIVecF <-
    dftExecute
      newPlan
      (DFTPlanID DFT1DG [numAngularFreq, numRadialFreq, cols, rows] [2, 3])
      filterPIVec
  return
    ( newPlan
    , fromUnboxed (Z :. numAngularFreq :. numRadialFreq :. cols :. rows) .
      VG.convert $
      filterVecF
    , fromUnboxed (Z :. numAngularFreq :. numRadialFreq :. cols :. rows) .
      VG.convert $
      filterPIVecF)

{-# INLINE frequencyDomainMultiply #-}
frequencyDomainMultiply ::
     (Source r (Complex Double))
  => Array U DIM4 (Complex Double)
  -> Array r DIM2 (Complex Double)
  -> Array D DIM4 (Complex Double)
frequencyDomainMultiply filterF input =
  R.traverse2
    filterF
    input
    const
    (\f1 f2 idx@(Z :. af :. rf :. i :. j) -> f1 idx * f2 (Z :. i :. j))

{-# INLINE convolvePinwheel #-}
convolvePinwheel ::
     (Source r (Complex Double))
  => DFTPlan
  -> Array U DIM4 (Complex Double)
  -> Array r DIM2 (Complex Double)
  -> IO (Array U DIM4 (Complex Double))
convolvePinwheel plan filterF input =
  fmap (fromUnboxed (extent filterF) . VG.convert) .
  dftExecute
    plan
    (DFTPlanID IDFT1DG (L.reverse . listOfShape . extent $ filterF) [2, 3]) .
  VG.convert . toUnboxed . computeUnboxedS $
  frequencyDomainMultiply filterF input

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE RecordWildCards #-}
import System.IO
import Control.Monad
import Data.Complex
import Data.List
import Data.ByteString (ByteString)

data Obj2D = Obj2D [String]
data Obj3D = Obj3D [String]

class Render2D a where
  render2d :: a -> [Complex Double]

class Render3D a where
  render3d :: a -> Obj3D

data Gear = Gear {
  gAlpha :: Double
, gTooth :: Double
, gModule :: Double
, gWidth :: Double
} deriving Show

data Bar = Bar {
  bWidth :: Double
} deriving Show

involute :: Double -> Double -> Complex Double
involute r theta = (r*(cos theta + theta * sin theta)) :+ (r*(sin theta - theta * cos theta))

yrev :: Complex Double -> Complex Double
yrev (a:+b) = (a:+(-b))

alpha :: Gear -> Double
alpha Gear{..} = 2 * pi * gAlpha / 360

radius :: Gear -> Double
radius g@Gear{..} = (cos (alpha g)) * radius' g

radius' :: Gear -> Double
radius' Gear{..} = gModule * gTooth / 2

instance Render2D Gear where
  render2d g@Gear{..} =
    let dat :: [Complex Double]
        dat = ((radius g) - 0.5) :+ 0 : map (\i-> involute (radius g) (i*2.0*pi/100.0))  [0..10]
        dat''' :: [Complex Double]
        dat''' = map (\i-> mkPolar (radius g) (i*2.0*pi/100.0))  [0..100]
        gear''  =  concat $
                   flip map [0..(gTooth-1)] $
                   \t -> map (\i -> i * (cis (t*2*pi/gTooth) )) $
                   dat ++
                   map (\t -> t * cis (0.9*pi/gTooth))
                   (reverse (map yrev dat))
    in gear''

instance Render3D Gear where
  render3d g@Gear{..} =
    let gear2d = render2d g
        gear2d' = intercalate "," $ map (\(r:+i) -> "[" ++ show r ++ "," ++ show i ++ "]") gear2d
    in Obj3D ["linear_extrude(height = " ++ show gWidth ++ ") { polygon( points= [",
              gear2d',
              "] );}"]

writePlot :: String -> [Complex Double] -> IO ()
writePlot file dat = do
  withFile file WriteMode $ \h -> do
    forM_ dat $ \(r:+i) -> do
      hPutStrLn h $ show r ++ " " ++ show i

write3DPlot :: String -> Obj3D -> IO ()
write3DPlot file (Obj3D dat) = do
  withFile file WriteMode $ \h -> do
    forM_ dat $ \str -> do
      hPutStr h $ str

main = do
  writePlot "involute.dat" $ render2d $ Gear 20 8 1 6
  write3DPlot "involute.scad" $ render3d $ Gear 20 8 1 6

max_stars_repo_name: junjihashimoto/raspberry-pi-box

max_issues_repo_name: junjihashimoto/raspberry-pi-box

max_forks_repo_name: junjihashimoto/raspberry-pi-box


text:
module Analyzer where

import Parser
import Control.Lens
import qualified Data.Map as M
import Numeric.LinearAlgebra

dictToMatrix :: M.Map User Int
             -> M.Map (User,User) Int
             -> (Matrix R,[User])
dictToMatrix ud uud = ((n><n) probs,map fst ucount)
  where
    n = length ud
    uucount = M.toList uud
    ucount = M.toList ud
    uutotal = [(u,u',n) | (u,n) <- ucount, (u',_) <- ucount]
    probs = combine uutotal uucount
    combine :: [(User,User,Int)] -> [((User,User),Int)] -> [R]
    combine _ [] = []
    combine [] as = map (const 0) as
    combine ((u,u',n):ts) olds@(((v,v'),k):news)
      | u == v && u' == v' = fromIntegral k/fromIntegral n : combine ts news
      | otherwise = 0 : combine ts olds

matrixP = uncurry dictToMatrix . (\(d,d') -> (fst <$> d,d')) <$> dictP

max_stars_repo_name: bnarnold/WAAnalyzer

max_issues_repo_name: bnarnold/WAAnalyzer

max_forks_repo_name: bnarnold/WAAnalyzer