text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# OPTIONS_GHC -Wall #-}

-- | Metric classes
module Plankton.Metric
  ( Signed(..)
  , Normed(..)
  , Metric(..)
  , Epsilon(..)
  , (≈)
  ) where

import Data.Complex (Complex(..))
import Plankton.Additive
import Plankton.Field
import Plankton.Multiplicative
import qualified Protolude as P
import Protolude
       (Bool(..), Double, Eq(..), Float, Int, Integer, Ord(..), ($), (&&))

-- | 'signum' from base is not an operator replicated in numhask, being such a very silly name, and preferred is the much more obvious 'sign'.  Compare with 'Norm' and 'Banach' where there is a change in codomain
--
-- > abs a * sign a == a
--
-- Generalising this class tends towards size and direction (abs is the size on the one-dim number line of a vector with its tail at zero, and sign is the direction, right?).
class (MultiplicativeUnital a) =>
      Signed a where
  sign :: a -> a
  abs :: a -> a

instance Signed Double where
  sign a =
    if a >= zero
      then one
      else negate one
  abs = P.abs

instance Signed Float where
  sign a =
    if a >= zero
      then one
      else negate one
  abs = P.abs

instance Signed Int where
  sign a =
    if a >= zero
      then one
      else negate one
  abs = P.abs

instance Signed Integer where
  sign a =
    if a >= zero
      then one
      else negate one
  abs = P.abs

-- | Like Signed, except the codomain can be different to the domain.
class Normed a b where
  size :: a -> b

instance Normed Double Double where
  size = P.abs

instance Normed Float Float where
  size = P.abs

instance Normed Int Int where
  size = P.abs

instance Normed Integer Integer where
  size = P.abs

instance (Multiplicative a, ExpField a, Normed a a) =>
         Normed (Complex a) a where
  size (rx :+ ix) = sqrt (rx * rx + ix * ix)

-- | distance between numbers
--
-- > distance a b >= zero
-- > distance a a == zero
-- > \a b c -> distance a c + distance b c - distance a b >= zero &&
-- >           distance a b + distance b c - distance a c >= zero &&
-- >           distance a b + distance a c - distance b c >= zero &&
class Metric a b where
  distance :: a -> a -> b

instance Metric Double Double where
  distance a b = abs (a - b)

instance Metric Float Float where
  distance a b = abs (a - b)

instance Metric Int Int where
  distance a b = abs (a - b)

instance Metric Integer Integer where
  distance a b = abs (a - b)

instance (Multiplicative a, ExpField a, Normed a a) =>
         Metric (Complex a) a where
  distance a b = size (a - b)

-- | todo: This should probably be split off into some sort of alternative Equality logic, but to what end?
class (AdditiveGroup a) =>
      Epsilon a where
  nearZero :: a -> Bool
  aboutEqual :: a -> a -> Bool
  positive :: (Eq a, Signed a) => a -> Bool
  positive a = a == abs a
  veryPositive :: (Eq a, Signed a) => a -> Bool
  veryPositive a = P.not (nearZero a) && positive a
  veryNegative :: (Eq a, Signed a) => a -> Bool
  veryNegative a = P.not (nearZero a P.|| positive a)

infixl 4 ≈

-- | todo: is utf perfectly acceptable these days?
(≈) :: (Epsilon a) => a -> a -> Bool
(≈) = aboutEqual

instance Epsilon Double where
  nearZero a = abs a <= (1e-12 :: Double)
  aboutEqual a b = nearZero $ a - b

instance Epsilon Float where
  nearZero a = abs a <= (1e-6 :: Float)
  aboutEqual a b = nearZero $ a - b

instance Epsilon Int where
  nearZero a = a == zero
  aboutEqual a b = nearZero $ a - b

instance Epsilon Integer where
  nearZero a = a == zero
  aboutEqual a b = nearZero $ a - b

instance (Epsilon a) => Epsilon (Complex a) where
  nearZero (rx :+ ix) = nearZero rx && nearZero ix
  aboutEqual a b = nearZero $ a - b

max_stars_repo_name: chessai/plankton

max_issues_repo_name: chessai/plankton

max_forks_repo_name: chessai/plankton


text:
-- Graphical and Network Models in Haskell.
-- Nicolas Kim
-- Department of Statistics, Carnegie Mellon University
-- Distributed under the MIT License.

import System.IO
import System.Random

import Adjacency
import Core
import Local

import Control.Monad
import Control.Monad.Random
import Control.Parallel (par, pseq)
-- import Control.Lens

import qualified Data.Map.Lazy as M
import qualified Data.List as L

import Numeric.LinearAlgebra (dispf)

-- import qualified Numeric.LinearAlgebra as LA

-- Some examples of graphs.
gr1 = Graph [1..3] [(1,2), (1,3), (2,3)] :: Graph Int
gr2 = Graph [1..2] [(1,2), (1,3), (2,3)] :: Graph Int
gr3 = completeGraph [1..10] :: Graph Int
gr4 = Graph [1..4] [(1,2), (3,4), (1,4)] :: Graph Int

kite = Graph [1..4] [(1,2), (1,3), (2,3), (2,4), (3,4)] :: Graph Int
wiki = Graph [1..6] [(1,2), (1,5), (2,5), (2,3), (3,4), (4,5), (4,6)] ::
  Graph Int
gr5 = Graph [1..3] [(1,2)] :: Graph Int
disconnected = Graph [1..4] [(1,2), (3,4)] :: Graph Int

-- Data structure for network models.
-- data Model = ERG { parameters :: [Double] } deriving (Show, Eq)

-- Makes an empty graph.
emptyGraph :: Graph a
emptyGraph = Graph [] []

-- Makes a complete simple graph from a list of nodes.
completeGraph :: (Ord a, Eq a) => [a] -> Graph a
completeGraph ns = Graph sns es
  where
    sns = L.sort ns
    es = completeEdges sns

-- List of all possible graphs on n nodes.
-- allGraphs :: (Ord a, Eq a) => [a] -> [Graph a] -> M.Map [Int] [Graph a]
-- allGraphs ns = L.foldl' () mapInit [1..n]
--   where
--     mapInit = M.fromList [(0, emptyGraph)]
--     n = (\x -> x*(x-1) `div` 2) $ length ns

-- Makes a valid graph object from a list of edges.
toGraph :: (Ord a, Eq a) => [(a, a)] -> Graph a
toGraph es = Graph (extractNodes es) es

-- Checks if the edges contain defined nodes.
validGraph :: (Eq a) => Graph a -> Bool
validGraph (Graph ns es) = subset (L.union (map fst es) (map snd es)) ns
  where
    subset xs ys = all (`elem` ys) xs

-- Removes redundant edges (e.g. (1,2) = (2,1)).
redundant :: (Ord a, Eq a) => [(a, a)] -> [(a, a)]
redundant es = L.intersect es (map f es)
  where
    f (x, y) = (min x y, max x y)

-- Makes a graph valid (see the `validGraph` function).
makeValid :: (Eq a) => Graph a -> Graph a
makeValid gr@(Graph ns es)
  | validGraph gr = gr
  | otherwise = Graph (L.union ens ns) es
    where
      ens = L.nub $ L.union (map fst es) (map snd es)

-- From a list of edges, give the sorted and pruned list of nodes
extractNodes :: (Ord a, Eq a) => [(a, a)] -> [a]
extractNodes es = L.nub $ L.sort allnodes
  where
    allnodes = (map fst es) ++ (map snd es)

-- The complete set of edges between nodes in a list. Nodes don't need to
-- be orderable. 22s
completeEdges' :: [a] -> [(a, a)]
completeEdges' ns = [ (fst x, fst y) | x <- zs, y <- zs, snd x < snd y ]
  where
    zs = zip ns [1..]

-- Add a set of nodes to a graph using their labels.
addNodes :: (Ord a) => Graph a -> [a] -> Graph a
addNodes (Graph ns es) nns = Graph (L.union ns nns) es

-- Add a set of edges to a graph.
addEdges :: (Ord a) => Graph a -> [(a, a)] -> Graph a
addEdges (Graph ns es) ees = Graph ns (L.union es ees)

-- -- Various graph statistics.
-- Computes the number of nodes.
nNode :: Graph a -> Int
nNode gr = length $ nodes gr

-- Computes the number of edges.
nEdge :: Graph a -> Int
nEdge gr = length $ edges gr

-- -- Monad-based generation of random nodes and edges, i.e. random graphs.
-- Erdos-Renyi Model
-- Include or don't include an edge?
pBool :: (RandomGen g) => Double -> Rand g Bool
pBool p = liftM (< p) $ getRandomR ((0, 1) :: (Double, Double))

-- Decide over a list of edges.
pBools :: (RandomGen g) => [Double] -> Rand g [Bool]
pBools ps = sequence $ map pBool ps

-- Given a list of Bools and any other list, return the list subsetted by
-- the list of Bools.
byBool :: [Bool] -> [a] -> [a]
byBool bs xs = [ snd zs | zs <- (zip bs xs), fst zs ]

-- Select elements from a list independently and with probabilities ps.
setEdges :: (RandomGen g) => [a] -> [Double] -> Rand g [a]
setEdges xs ps = liftM (flip byBool xs) $ pBools ps

-- Generate an Erdos-Renyi random graph; nodes ns, and edge probabilities
-- ps.
erdosGen :: (Ord a, RandomGen g) => [a] -> [Double] -> Rand g (Graph a)
erdosGen ns ps = liftM (Graph ns) $ setEdges (completeEdges ns) ps

-- Graphon Model: generates a w-random graph.
graphonGen :: (Ord a, RandomGen g) => [a] -> (Double -> Double -> Double)
  -> Rand g (Graph a)
graphonGen ns w = liftM (Graph ns) es
  where
    us = liftM (take (length ns)) $ getRandomRs ((0, 1) :: (Double, Double))
    es = (liftM (applyUpper w) us) >>= (setEdges (completeEdges ns))

-- Graphon Model: generates a w-random graph. Parallel.
graphonGen' :: (Ord a, RandomGen g) => [a] -> (Double -> Double -> Double)
  -> Rand g (Graph a)
graphonGen' ns w = liftM2 pseq (liftM2 par (liftM force es) (liftM force
                               es')) (liftM (Graph ns) (liftM2 (++) es
                               es'))
  where
    us = liftM (take (length ns)) $ getRandomRs ((0, 1) :: (Double, Double))
    cs = completeEdges ns
    n = length ns
    cshalf = splitAt (((n * (n-1))+2) `div` 4) cs
    es = (liftM (applyUpper w) us) >>= (setEdges (fst cshalf))
    es' = (liftM (applyUpper w) us) >>= (setEdges (snd cshalf))

-- From "Real World Haskell". Used in the parallel version of the code for
-- generating w-random graphs, `graphonGen'`.
force :: [a] -> ()
force xs = go xs `pseq` ()
  where go (_:xs) = go xs
        go [] = 1

-- Apply a function to the upper triangle of an array.
applyUpper :: (a -> a -> b) -> [a] -> [b]
applyUpper f xs = [ f (fst x) (fst y) | x <- zl, y <- zl, snd x < snd y ]
  where
    zl = zip xs ([1..] :: [Int])

sblock :: Double -> Double -> Double
sblock x y
  | x < 0.5 && y < 0.5 = 0.8
  | x < 0.5 || y < 0.5 = 0.2
  | otherwise          = 0.8

-- Need to create a show instance for this (or rewrite it).
round' :: (RealFrac a, RealFrac b, Integral b) => Int -> a -> b
round' x = (/ 10^x) . round . (* 10^x)

main :: IO ()
main = do
    -- values <- evalRandIO $ graphonGen [1..10] (\x y -> (x+y)/2)
    let n = 100 :: Int
    gr <- evalRandIO $ graphonGen [1..n] sblock
    putStrLn . show . length $ community (toMap gr) 1
    -- let spannum = spanTreeCount values
    -- putStrLn . (\x -> "Number of spanning trees: " ++ x) $ show spannum
    -- putStrLn . show $ spectralCluster' values
    -- putStrLn . (\x -> "Guess for optimal lambda: " ++ x) . show
    --  $ lambdaSelect values
    -- putStrLn . (\x -> "This is " ++ x ++ "% of the true value." ) . show
    --   . round' 4 . ((/) $ fromIntegral n^(n-2)) $ fromIntegral spannum

ext: hs

max_stars_repo_name: kimolas/causal-haskell

max_issues_repo_name: kimolas/causal-haskell

max_forks_repo_name: kimolas/causal-haskell


text:
-- | An ad-hoc representation of matrices that saves allocations
-- and speeds up computing gradients. The major improvement comes
-- from avoiding, often, the construction of some matrices,
-- e.g., outer products or repeated columns that naturally occur
-- in delta-expressions. Minor improvements come from better fusion
-- (via vector or hmatrix/blas/lapack) and so less allocation
-- of intermediate matrix results of arithmetic and other operations.
--
-- The latter improvements are very likely tied to the vagaries
-- of hmatrix (or, less likely, vector or blas/lapack) that work underneath
-- and apparently conspire to fuse some matrix operations but not others.
-- That would explain why vector-based MNIST nn allocates not much less
-- than matrix-based MNIST, despite producing much larger delta-expressions.
-- Matrix operations probably fuse worse, even though they are really
-- vector operations on the underlying vector representation,
-- because the dimensions book-keeping of the matrix comes in the way
-- and also because some operations work on columns, against the grain
-- of the representation.
--
-- For fully connected MNIST nns 10 times larger than usual, vector-based
-- implementation allocates similarly and runs faster than matrix-based,
-- despite the latter spending 8 times less in GC. probably because
-- matrices don't fit in cache at this size and vectors still do.
module HordeAd.Internal.MatrixOuter
  ( MatrixOuter (..)
  , emptyMatrixOuter, nullMatrixOuter, rows, cols
  , convertMatrixOuter, convertMatrixOuterOrNull
  , toRows, toColumns, plus, multiplyWithOuter
  , konst, asRow, asColumn, rowAppend, columnAppend
  , takeRows, takeColumns, dropRows, dropColumns
  , transpose, matMul, flipud, fliprl
  ) where

import Prelude

import           Control.Exception.Assert.Sugar
import           Data.Maybe (fromJust)
import qualified Data.Vector.Generic as V
import           Numeric.LinearAlgebra (Matrix, Numeric, Vector)
import qualified Numeric.LinearAlgebra as HM

-- | A representation of a matrix as a product of a basic matrix
-- and an outer product of two vectors. Each component defaults to ones.
data MatrixOuter r = MatrixOuter (Maybe (Matrix r))
                                 (Maybe (Vector r)) (Maybe (Vector r))

emptyMatrixOuter :: MatrixOuter r
emptyMatrixOuter = MatrixOuter Nothing Nothing Nothing

nullMatrixOuter :: MatrixOuter r -> Bool
nullMatrixOuter (MatrixOuter Nothing Nothing Nothing) = True
nullMatrixOuter _ = False

rows :: Numeric r => MatrixOuter r -> Int
rows (MatrixOuter (Just m) _ _) = HM.rows m
rows (MatrixOuter Nothing (Just c) _) = V.length c
rows _ = error "rows: dimensions can't be determined"

cols :: Numeric r => MatrixOuter r -> Int
cols (MatrixOuter (Just m) _ _) = HM.cols m
cols (MatrixOuter Nothing _ (Just r)) = V.length r
cols _ = error "cols: dimensions can't be determined"

convertMatrixOuter :: (Numeric r, Num (Vector r)) => MatrixOuter r -> Matrix r
convertMatrixOuter (MatrixOuter (Just m) Nothing Nothing) = m
convertMatrixOuter (MatrixOuter (Just m) (Just c) Nothing) = m * HM.asColumn c
  -- This @asColumn@ is safe, because it's expanded via @m@.
  --
  -- Beware, depending on context, @m * outer c (konst 1 (cols m))@
  -- may allocate much less and perhaps @fromRows . toRowsMatrixOuter@
  -- may be even better; that's probably vector being picky about fusing.
  -- It doesn't matter if @m@ comes first.
convertMatrixOuter (MatrixOuter (Just m) Nothing (Just r)) = m * HM.asRow r
  -- This @asRow@ is safe, because it's expanded via @m@.
convertMatrixOuter (MatrixOuter (Just m) (Just c) (Just r)) = m * HM.outer c r
convertMatrixOuter (MatrixOuter Nothing (Just c) (Just r)) = HM.outer c r
convertMatrixOuter _ =
  error "convertMatrixOuter: dimensions can't be determined"

convertMatrixOuterOrNull :: (Numeric r, Num (Vector r))
                         => MatrixOuter r -> Matrix r
convertMatrixOuterOrNull (MatrixOuter Nothing Nothing Nothing) = HM.fromRows []
convertMatrixOuterOrNull m = convertMatrixOuter m

toRows :: (Numeric r, Num (Vector r)) => MatrixOuter r -> [Vector r]
toRows (MatrixOuter (Just m) Nothing Nothing) = HM.toRows m
toRows (MatrixOuter (Just m) mc Nothing) =
  maybe id (zipWith HM.scale . V.toList) mc $ HM.toRows m
toRows (MatrixOuter (Just m) mc (Just r)) =
  maybe (map (r *)) (zipWith (\s row -> r * HM.scale s row) . V.toList) mc
  $ HM.toRows m
toRows (MatrixOuter Nothing (Just c) (Just r)) =
  map (`HM.scale` r) $ V.toList c
toRows _ = error "toRows: dimensions can't be determined"

toColumns :: (Numeric r, Num (Vector r)) => MatrixOuter r -> [Vector r]
toColumns (MatrixOuter (Just m) Nothing Nothing) = HM.toColumns m
toColumns (MatrixOuter (Just m) Nothing mr) =
  maybe id (zipWith HM.scale . V.toList) mr $ HM.toColumns m
toColumns (MatrixOuter (Just m) (Just c) mr) =
-- traceShow (sizeMatrixOuter (MatrixOuter (Just m) (Just c) mr)) $
  maybe (map (c *)) (zipWith (\s col -> c * HM.scale s col) . V.toList) mr
  $ HM.toColumns m
toColumns (MatrixOuter Nothing (Just c) (Just r)) =
  map (`HM.scale` c) $ V.toList r
toColumns _ = error "toColumns: dimensions can't be determined"

size :: Numeric r => MatrixOuter r -> (Maybe (Int, Int), Maybe Int, Maybe Int)
size (MatrixOuter mm mc mr) = (HM.size <$> mm, HM.size <$> mc, HM.size <$> mr)

plus :: (Numeric r, Num (Vector r))
     => MatrixOuter r -> MatrixOuter r -> MatrixOuter r
plus o1 o2 =
  let !o = convertMatrixOuter o1 + convertMatrixOuter o2
  in MatrixOuter (Just o) Nothing Nothing
    -- TODO: Here we allocate up to 5 matrices, but we should allocate one
    -- and in-place add to it and multiply it, etc., ideally using raw FFI
    -- or at least vector streams/bundles; an easy option to have only
    -- one allocation is to do this pointwise, but then the speedup, if any,
    -- from blas/lapack does not apply

multiplyWithOuter :: (Numeric r, Num (Vector r))
                  => Matrix r -> MatrixOuter r -> MatrixOuter r
multiplyWithOuter k o@(MatrixOuter mm mc mr) =
  assert (HM.size k
          == maybe (HM.size $ fromJust mc, HM.size $ fromJust mr) HM.size mm
         `blame` (HM.size k, size o)) $
  let !m = maybe k (k *) mm
  in MatrixOuter (Just m) mc mr

konst :: Numeric r => r -> Int -> Int -> MatrixOuter r
konst x nrows ncols =
  let !c = HM.konst x nrows
      !r = HM.konst x ncols
  in MatrixOuter Nothing (Just c) (Just r)

asRow :: Numeric r => Vector r -> Int -> MatrixOuter r
asRow !r nrows =
  let !c = HM.konst 1 nrows
  in MatrixOuter Nothing (Just c) (Just r)

asColumn :: Numeric r => Vector r -> Int -> MatrixOuter r
asColumn !c ncols =
  let !r = HM.konst 1 ncols
  in MatrixOuter Nothing (Just c) (Just r)

rowAppend :: (Numeric r, Num (Vector r))
          => MatrixOuter r -> MatrixOuter r -> MatrixOuter r
rowAppend o1 o2 =
  let !o = convertMatrixOuter o1 HM.=== convertMatrixOuter o2
  in MatrixOuter (Just o) Nothing Nothing

columnAppend :: (Numeric r, Num (Vector r))
             => MatrixOuter r -> MatrixOuter r -> MatrixOuter r
columnAppend o1 o2 =
  let !o = convertMatrixOuter o1 HM.||| convertMatrixOuter o2
  in MatrixOuter (Just o) Nothing Nothing

takeRows :: Numeric r => Int -> MatrixOuter r -> MatrixOuter r
takeRows k (MatrixOuter mm mc mr) =
  MatrixOuter (HM.takeRows k <$> mm)
              (V.take k <$> mc)
              mr

takeColumns :: Numeric r => Int -> MatrixOuter r -> MatrixOuter r
takeColumns k (MatrixOuter mm mc mr) =
  MatrixOuter (HM.takeColumns k <$> mm)
              mc
              (V.take k <$> mr)

dropRows :: Numeric r => Int -> MatrixOuter r -> MatrixOuter r
dropRows k (MatrixOuter mm mc mr) =
  MatrixOuter (HM.dropRows k <$> mm)
              (V.drop k <$> mc)
              mr

dropColumns :: Numeric r => Int -> MatrixOuter r -> MatrixOuter r
dropColumns k (MatrixOuter mm mc mr) =
  MatrixOuter (HM.dropColumns k <$> mm)
              mc
              (V.drop k <$> mr)

transpose :: Numeric r => MatrixOuter r -> MatrixOuter r
transpose (MatrixOuter mm mc mr) = MatrixOuter (HM.tr' <$> mm) mr mc

matMul :: (Numeric r, Num (Vector r))
       => MatrixOuter r -> MatrixOuter r -> MatrixOuter r
matMul m1 m2 =
  MatrixOuter (Just $ convertMatrixOuter m1 HM.<> convertMatrixOuter m2)
              Nothing Nothing

flipud :: Numeric r => MatrixOuter r -> MatrixOuter r
flipud (MatrixOuter mm mc mr) =
  MatrixOuter (HM.flipud <$> mm) mr (V.reverse <$> mc)

fliprl :: Numeric r => MatrixOuter r -> MatrixOuter r
fliprl (MatrixOuter mm mc mr) =
  MatrixOuter (HM.fliprl <$> mm) (V.reverse <$> mr) mc

max_stars_repo_name: Mikolaj/horde-ad

max_issues_repo_name: Mikolaj/horde-ad

max_forks_repo_name: Mikolaj/horde-ad


text:
{-# LANGUAGE NoMonomorphismRestriction #-}
module Widget.Layout
  ( Layout(..)
  , mkLayout
  , arrangeLayout
  , resetLayout
  ) where

import Import
import File
import Config

import Statistics.Distribution.Empirical (levelEmpiricalDistribution, cumulativeEmpiricalDistribution)

import Widget.Title
import Widget.Splits
import Widget.Times
import Widget.Total
import Widget.Percentile

data Layout = L Splits Percentile Times Title Total

instance Widget Layout where
  update = updateLayout
  handle = handleLayout
  redraw = redrawLayout
  window = error "Called `window` on Layout"
  move   = error "Called `move` on Layout"
  resize = error "Called 'resize' on Layout"

mkLayout :: Config -> FileFormat -> Curses Layout
mkLayout (Config title levelNames) file = do
  -- Window setup
  let levelBestTimes  = map snd $ levelBests file
      cumulativeTimes = map snd $ cumulativeEmpiricalDistribution file
      levelDistrs     = map snd $ levelEmpiricalDistribution file
      percentileData  = zip levelDistrs levelBestTimes

  layout <- L <$> (mkSplits levelNames          <$> newWindow')
              <*> (mkPercentile percentileData  <$> newWindow')
              <*> (mkTimes file                 <$> newWindow')
              <*> (mkTitle file title           <$> defaultWindow)
              <*> (mkTotal file cumulativeTimes <$> newWindow')

  arrangeLayout layout

  return layout

updateLayout :: NominalDiffTime -> Layout -> Curses Layout
updateLayout d (L s p times title to) =
  let f = update d
  in L <$> f s
       <*> f p
       <*> f times
       <*> f title
       <*> f to

handleLayout :: TimerAction -> Layout -> Curses Layout
handleLayout ta (L s p times title to) =
  let f = handle ta
  in L <$> f s
       <*> f p
       <*> f times
       <*> f title
       <*> f to

redrawLayout :: Layout -> Curses ()
redrawLayout (L s p times title to) = do
  redraw s
  redraw p
  redraw times
  redraw title
  redraw to

newWindow' :: Curses Window
newWindow' = newWindow 1 1 0 0

arrangeLayout :: Layout -> Curses ()
arrangeLayout (L splits percentiles times title total) = do
  (rows, columns) <- screenSize

  let titleRows  = 2
      totalsRows = 4
      splitsRows = rows - (titleRows + totalsRows)

      timesColumns      = 10
      percentileColumns = 5
      splitsColumns     = columns - (timesColumns + percentileColumns)

  move   title 0 0
  resize title titleRows columns
  redraw title

  move   splits titleRows 0
  resize splits splitsRows splitsColumns
  redraw splits

  move   percentiles titleRows splitsColumns
  resize percentiles splitsRows percentileColumns
  redraw percentiles

  move   times titleRows (splitsColumns + percentileColumns)
  resize times splitsRows timesColumns
  redraw times

  move   total (titleRows + splitsRows) 0
  resize total totalsRows columns
  redraw total

resetLayout :: Config -> FileFormat -> Layout -> Layout
resetLayout (Config title levelNames) file (L splits percentiles times titleW total) =
  let levelBestTimes  = map snd $ levelBests file
      cumulativeTimes = map snd $ cumulativeEmpiricalDistribution file
      levelDistrs     = map snd $ levelEmpiricalDistribution file
      percentileData  = zip levelDistrs levelBestTimes
  in L (mkSplits levelNames $ window splits)
       (mkPercentile percentileData $ window percentiles)
       (mkTimes file $ window times)
       (mkTitle file title $ window titleW)
       (mkTotal file cumulativeTimes $ window total)

max_stars_repo_name: argiopetech/timer

max_issues_repo_name: argiopetech/timer

max_forks_repo_name: argiopetech/timer


text:
module Grid
    ( Grid
    , standardInterval
    , makeInterval
    , startOfInterval
    , endOfInterval
    , elementSize
    , nElements
    ) where

import Numeric.LinearAlgebra (R)

data Grid = Grid R R R Int
  deriving (Eq, Show)

standardInterval :: Grid
standardInterval = makeInterval 0.0 1.0 10

-- | Specify a 1D Intervall [start, end]
makeInterval :: R -> R -> Int -> Grid
makeInterval start end n
  | start < end && n > 0 = Grid start end ((end - start) / fromIntegral n) n
  | otherwise = undefined

startOfInterval :: Grid -> R
startOfInterval (Grid start _ _ _) = start

endOfInterval :: Grid -> R
endOfInterval (Grid _ end _ _) = end

elementSize :: Grid -> R
elementSize (Grid _ _ h _) = h

nElements :: Grid -> Int
nElements (Grid _ _ _ n) = n

max_stars_repo_name: bschwb/haskfem

max_issues_repo_name: bschwb/haskfem

max_forks_repo_name: bschwb/haskfem


text:
{- |

Module      :  Numeric.GSL
Copyright   :  (c) Alberto Ruiz 2006-7
License     :  GPL-style

Maintainer  :  Alberto Ruiz (aruiz at um dot es)
Stability   :  provisional
Portability :  uses -fffi and -fglasgow-exts

This module reexports all available GSL functions.

The GSL special functions are in the separate package hmatrix-special.

-}

module Numeric.GSL (
  module Numeric.GSL.Integration
, module Numeric.GSL.Differentiation
, module Numeric.GSL.Fourier
, module Numeric.GSL.Polynomials
, module Numeric.GSL.Minimization
, module Numeric.GSL.Root
, module Numeric.GSL.ODE
, module Numeric.GSL.Fitting
, module Data.Complex
, setErrorHandlerOff
) where

import Numeric.GSL.Integration
import Numeric.GSL.Differentiation
import Numeric.GSL.Fourier
import Numeric.GSL.Polynomials
import Numeric.GSL.Minimization
import Numeric.GSL.Root
import Numeric.GSL.ODE
import Numeric.GSL.Fitting
import Data.Complex

-- | This action removes the GSL default error handler (which aborts the program), so that
-- GSL errors can be handled by Haskell (using Control.Exception) and ghci doesn't abort.
foreign import ccall unsafe "GSL/gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO ()

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
module LayersSpec
  ( spec
  ) where

import           Control.Monad
import qualified Data.Set                as Set
import           Data.Tuple.Extra
import           Debug.Trace
import           Foundation.Monad
import           Layers
import           Learn
import           Mnist
import           Numeric.LinearAlgebra   as A
import           Synapse
import           Test.Hspec
import           Test.Hspec.QuickCheck   (prop)
import           Test.QuickCheck
import           Test.QuickCheck.Gen
import           Test.QuickCheck.Monadic

runTest :: IO ()
runTest = hspec spec

spec :: Spec
spec = do
  describe "predict" propsPredict
  describe "evaluate" propsEvaluate
  describe "Affine" propsAffine
  describe "Sigmoid" propsSigmoid
  describe "ReLU" propsReLU
  describe "Joined" propsJoined
  describe "Matrix" propsMatrix

propsPredict =
  prop "index" $
  forAll (genVectorN 100) $ \v ->
    let r = predict ReLUForward v
        a = maxIndex v
     in r `shouldBe` a

propsEvaluate = do
  prop "100%" $
    forAll (vectorOf 100 $ genVectorN 100) $ \vs ->
      let ns = map maxIndex vs
          r = evaluate ReLUForward $ zip ns vs
       in r `shouldBe` 1
  prop "0%" $
    forAll (vectorOf 100 $ genVectorN 100) $ \vs ->
      let ns = map ((+ 1) . maxIndex) vs
          r = evaluate ReLUForward $ zip ns vs
       in r `shouldBe` 0
  prop "50%" $
    forAll (vectorOf 100 $ genVectorN 100) $ \vs ->
      let (vA, vB) = splitAt 50 vs
          nA = map maxIndex vA
          nB = map ((+ 1) . maxIndex) vB
          ns = nA ++ nB
          r = evaluate ReLUForward $ zip ns vs
       in r `shouldBe` 0.5

propsAffine = do
  prop "forward" $
    forAll
      (do a <- choose (10, 100)
          b <- choose (10, 100)
          c <- choose (10, 100)
          f <- genAffine a b
          s <- genMatrix c a
          return (f, s)) $ \(af, x) ->
      let AffineForward weights bias = af
          (AffineBackward rW rB rX, r) = forward af x
          e = affinem weights bias x
       in (rW, rB, rX, r) `shouldBe` (weights, bias, x, e)
  prop "backward" $
    forAll
      (do a <- choose (10, 100)
          b <- choose (10, 100)
          c <- choose (10, 100)
          d <- genMatrix c b
          f <- genAffine a b
          s <- genMatrix c a
          return (f, s, d)) $ \(af, x, d) ->
      let AffineForward w b = af
          rate = 1.0
          (AffineForward w' b', d') = backward rate (AffineBackward w b x) d
          (dx, dw, db) = affinemBackward w x d
          rW = w - dw
          rB = b - db
       in (w', b', d') `shouldBe` (rW, rB, dx)
  prop "backward with rate" $
    forAll
      (do a <- choose (10, 100)
          b <- choose (10, 100)
          c <- choose (10, 100)
          r <- arbitrary `suchThat` (\a -> 0 < a && a < 1)
          d <- genMatrix c b
          f <- genAffine a b
          s <- genMatrix c a
          return (r, f, s, d)) $ \(rate, af, x, d) ->
      let AffineForward w b = af
          (AffineForward w' b', d') = backward rate (AffineBackward w b x) d
          (dx, dw, db) = affinemBackward w x d
          rW = w - A.scale rate dw
          rB = b - A.scale rate db
       in (w', b', d') `shouldBe` (rW, rB, dx)

propsSigmoid = do
  prop "forward" $
    forAll
      (do a <- choose (10, 100)
          b <- choose (10, 100)
          genMatrix a b) $ \x ->
      let (SigmoidBackward y, r) = forward SigmoidForward x
          e = sigmoidm x
       in (y, r) `shouldBe` (r, e)
  prop "backward" $
    forAll
      (do a <- choose (10, 100)
          b <- choose (10, 100)
          d <- genMatrix a b
          y <- genMatrix a b
          return (y, d)) $ \(y, d) ->
      let (SigmoidForward, d') = backward 0 (SigmoidBackward y) d
          e = sigmoidBackward y d
       in d' `shouldBe` e

propsReLU = do
  prop "forward" $
    forAll
      (do a <- choose (10, 100)
          b <- choose (10, 100)
          genMatrix a b) $ \x ->
      let (ReLUBackward y, r) = forward ReLUForward x
          e = relum x
       in (y, r) `shouldBe` (x, e)
  prop "backward" $
    forAll
      (do a <- choose (10, 100)
          b <- choose (10, 100)
          x <- genMatrix a b
          d <- genMatrix a b
          return (x, d)) $ \(x, d) ->
      let (ReLUForward, d') = backward 0 (ReLUBackward x) d
          e = relumBackward x d
       in d' `shouldBe` e

propsJoined = do
  prop "forward" $
    forAll
      (do a <- choose (10, 100)
          b <- choose (10, 100)
          c <- choose (10, 100)
          f <- genAffine a b
          x <- genMatrix c a
          return (f, x)) $ \(fL1, x) ->
      let fL2 = SigmoidForward
          fJ = fL1 ~> fL2
          (JoinedBackwardLayer bL1 bL2, y) = forward fJ x
          (bL1', yL1) = forward fL1 x
          (bL2', yL2) = forward fL2 yL1
       in (bL1, bL2, y) `shouldBe` (bL1', bL2', yL2)
  prop "backward" $
    forAll
      (do a <- choose (10, 100)
          b <- choose (10, 100)
          c <- choose (10, 100)
          d <- genMatrix c b
          y <- genMatrix c b
          x <- genMatrix c a
          fL1 <- genAffine a b
          let AffineForward weights bias = fL1
          let bL1 = AffineBackward weights bias x
          return (bL1, SigmoidBackward y, d)) $ \(bL1, bL2, d) ->
      let bJ = bL1 <~ bL2
          (JoinedForwardLayer fL1 SigmoidForward, bD) = backward 1.0 bJ d
          (SigmoidForward, dL2) = backward 1.0 bL2 d
          (fL1', dL1) = backward 1.0 bL1 dL2
       in (fL1, bD) `shouldBe` (fL1', dL1)

propsMatrix = do
  prop "scale" $
    forAll genElements $ \(a, b, x, y) ->
      let m = b `A.scale` toMatrix x y a
          m' = toMatrix x y (a * b)
       in m `shouldBe` m'
  prop "multiply" $
    forAll genElements $ \(a, b, x, y) ->
      let m = toMatrix x y b * toMatrix x y a
          m' = toMatrix x y (a * b)
       in m `shouldBe` m'
  prop "subtract" $
    forAll genElements $ \(a, b, x, y) ->
      let m = toMatrix x y a - toMatrix x y b
          m' = toMatrix x y (a - b)
       in m `shouldBe` m'
  prop "add" $
    forAll genElements $ \(a, b, x, y) ->
      let m = toMatrix x y a + toMatrix x y b
          m' = toMatrix x y (a + b)
       in m `shouldBe` m'
  where
    toMatrix x y a = fromLists $ replicate x $ replicate y a
    genElements = do
      a <- arbitrary
      b <- arbitrary
      x <- choose (10, 100)
      y <- choose (10, 100)
      return (a :: R, b, x, y)

genVectorN :: Int -> Gen (Vector R)
genVectorN n = fmap fromList $ vectorOf n $ choose (1, 100)

genAffine :: Int -> Int -> Gen (ForwardLayer R)
genAffine nIn nOut = do
  weights <- fmap fromLists $ vectorOf nIn $ vectorOf nOut arbitrary
  bias <- fromList <$> vectorOf nOut arbitrary
  return $ AffineForward weights bias

genMatrix :: Int -> Int -> Gen (Matrix R)
genMatrix nRows nCols =
  fmap fromLists $ vectorOf nRows $ vectorOf nCols arbitrary


text:
module Utils where

import           Control.Monad         (ap, liftM, liftM2)
import           Data.Char             (isNumber)
import           Data.List             (genericLength)
import qualified Data.Vector           as V
import qualified Data.Vector.Unboxed   as U
import           Numeric.LinearAlgebra (Matrix, cmap, cols, dropColumns,
                                        flatten, fromRows, loadMatrix,
                                        takeColumns, toList, toRows)
import qualified Numeric.LinearAlgebra as L
import           System.FilePath       (splitFileName)
import           System.Posix          (fileExist)

usage :: String
usage = "Usage: ./init_present <file> <Graph nbd size> <PCA nbd size>"

errorMsg1 :: String
errorMsg1 = "You did not enter two positive integers.\n"

errorMsg2 :: String
errorMsg2 = "You did not enter two positive integer arguments.\n"

graphMsg1 :: String
graphMsg1 = "The graph is not fully connected. Please enter larger integers!"

fileMsg :: String
fileMsg = "You have not provided an existent file. Please try again."

getFile :: FilePath -> Bool -> FilePath
getFile file pCond = if pCond
                       then file
                       else error fileMsg

fileArgs :: [String] -> IO Bool
fileArgs = (=<<) fileExist . return . head

intCheck :: Foldable t => t Char -> Bool
intCheck x = (not . null $ x) && foldr ((&&) . isNumber) True x

gtZero :: [String] -> Bool
gtZero = all ((> 0) . (\x -> read x :: Int))

procMNIST :: [Double] -> Matrix Double -> Matrix Double
procMNIST selection mat = images
  where
    images = fromRows [x | (x, c) <- zip xs cs
                         , or $ sequenceA eqfns c]
    xs = toRows . cmap (\x -> 255 - x :: Double) . mSel takeColumns $ mat
    cs = toList $ flatten . mSel dropColumns $ mat
    eqfns = [(== x) | x <- selection]
    mSel f x = f ((+ (-1)) $ cols x) x

checkArgs :: [String] -> (Int,Int)
checkArgs args = if ((== (3 :: Integer)) . genericLength) args
                   then let val = tail args
                        in if all intCheck val && gtZero val
                             then head . ap zip tail . Prelude.map read $ val -- yeah...
                             else error (errorMsg1 ++ usage)
                   else error (errorMsg2 ++ usage)

checkFile :: [FilePath] -> IO FilePath
checkFile args = if not . null $ args
                   then liftM2 getFile <$> return . head <*> fileArgs $ args
                   else error usage

getReqMatrix :: FilePath -> IO (Matrix Double)
getReqMatrix file = if (snd . splitFileName $ file) == "mnist.txt"
                      then liftM (procMNIST [4]) $ loadMatrix file
                      else loadMatrix file

convertDataset :: Matrix Double -> V.Vector (U.Vector Double)
convertDataset = V.map (U.fromList . L.toList) . V.fromList . L.toRows

rowSize :: Matrix Double -> Int
rowSize = L.rows

max_stars_repo_name: emmanueldenloye/manifoldRNC

max_issues_repo_name: emmanueldenloye/manifoldRNC

max_forks_repo_name: emmanueldenloye/manifoldRNC


text:
{-# LANGUAGE BinaryLiterals #-}
{-# LANGUAGE OverloadedLists #-}

module Day03 where

import Data.Bits (Bits (complement, setBit, testBit, (.&.)))
import qualified Data.Text as T
import qualified Data.Vector.Storable as Vec
import Motif (count)
import Numeric.LinearAlgebra.Data (I, Vector)
import qualified Relude.Unsafe as Unsafe
import Utils (binaryLines, showSolutions, elimination)

-- least significant bit first
-- >>> bitsAsInts 0b101
-- [1,0,1,0,0]
bitsAsInts :: Int -> Word32 -> Vector I
bitsAsInts bitLength i =
  Vec.generate
    bitLength
    ( \idx ->
        if i `testBit` idx then 1 else 0
    )

-- count the number of bits in each place value of a list of binary numbers
bitCounts :: Int -> [Word32] -> Vector I
bitCounts bitLength = sum . map (bitsAsInts bitLength)

-- >>> boolVecToBitsLE [False, True, True]
-- 6
boolVecToBitsLE :: Vector Bool -> Word32
boolVecToBitsLE = Vec.ifoldl' step 0
  where
    step :: Word32 -> Int -> Bool -> Word32
    step acc idx b
      | b = acc `setBit` idx
      | otherwise = acc

calculateGamma :: Int -> [Word32] -> Word32
calculateGamma bitLength is =
  boolVecToBitsLE . Vec.map isMostCommonBit . bitCounts bitLength $ is
  where
    len = fromIntegral $ length is

    isMostCommonBit :: I -> Bool
    isMostCommonBit cnt = cnt > (len `div` 2)

exampleReport :: [Word32]
exampleReport =
  [ 0b00100,
    0b11110,
    0b10110,
    0b10111,
    0b10101,
    0b01111,
    0b00111,
    0b11100,
    0b10000,
    0b11001,
    0b00010,
    0b01010
  ]

data RatingType = Oxygen | CO2

-- >>> calculateRating 5 Oxygen exampleReport
-- 23
-- >>> calculateRating 5 CO2 exampleReport
-- 10

calculateRating :: Int -> RatingType -> [Word32] -> Word32
calculateRating bitLength ratingType is = final
  where
    Just final =
      evalState
        (elimination (elimByBitCriteria ratingType) is)
        (bitLength - 1)

elimByBitCriteria :: RatingType -> [Word32] -> State Int [Word32]
elimByBitCriteria ratingType candidates = do
  bitIndex <- get
  let ones = count (`testBit` bitIndex) candidates
      zeros = length candidates - ones
      theBit = case ratingType of
        Oxygen -> ones >= zeros
        CO2 -> ones < zeros
      remaining = filter ((== theBit) . (`testBit` bitIndex)) candidates
  put $ bitIndex - 1
  pure remaining

solve :: Text -> Text
solve input = showSolutions p1 p2
  where
    Just is = binaryLines input

    bitLength :: Int
    bitLength = T.length . Unsafe.head . T.lines $ input

    gamma = calculateGamma bitLength is
    -- we mask to ensure irrelevant high bits are not set by complement
    nOnes n = 2 ^ n - 1
    mask = nOnes bitLength
    epsilon = complement gamma .&. mask
    p1 = gamma * epsilon

    oxygenRating = calculateRating bitLength Oxygen is
    co2Rating = calculateRating bitLength CO2 is
    p2 = oxygenRating * co2Rating


text:
module Main (main) where

import Neural.Matrix
import Data.Ord (comparing)
import Data.List (scanl, foldl', maximumBy)
import Numeric.LinearAlgebra (toList, vector)
import Control.Monad (forM_)
import Codec.Compression.GZip (decompress)
import qualified Data.ByteString.Lazy as BS

getImage s n = fromIntegral . BS.index s . (n*28^2 + 16 +) <$> [0..28^2 - 1]
getX     s n = (vector (getImage s n)) / 256
getLabel s n = fromIntegral $ BS.index s (n + 8)
getY     s n = vector $ fromIntegral . fromEnum . (getLabel s n ==) <$> [0..9]

printImage imgs = putStrLn . unlines . take 28 . map (take 28) . iterate (drop 28)
                           . map pixel2Char . getImage imgs
  where pixel2Char n = let s = " ·:o◍O0@" in s !! (fromIntegral n * length s `div` 256)

display d prob = show d ++ ": " ++ show prob

bestOf :: [Float] -> Int
bestOf = fst . maximumBy (comparing snd) . zip [0..]

small = [[0..99], [100..299], [300..599], [600..999]]
medium = [[0..299], [300..999], [1000..2199], [2200..4599]]
large = [[0..999], [1000..2999], [3000..5999], [6000..9999]]
rapid = [[0..299], [300..599], [600..899], [900..1199], [1200..1499], [1500..1799], [1800..2099]]

main :: IO ()
main = do
  [trainI, trainL, testI, testL] <- mapM ((decompress  <$>) . BS.readFile . ("examples/mnistData/"++))
                                    ["train-images-idx3-ubyte.gz"
                                    ,"train-labels-idx1-ubyte.gz"
                                    ,"t10k-images-idx3-ubyte.gz"
                                    ,"t10k-labels-idx1-ubyte.gz"]
  net <- initNet 1 [784, 30, 10] [relu, relu]
  let n = 42
  printImage testI n

  let
    example = getX testI n
    randomUpdate net x = updateNet categoricalCrossEntropy 0.002 (getX trainI x) (getY trainL x) net
    trainingNet = scanl (foldl' randomUpdate) net large
    trainedNet = last trainingNet
  forM_ trainingNet $ putStrLn . unlines . zipWith display [0..9] . softmax
                               . toList . snd . last . forwardPass example

max_stars_repo_name: Ultramann/hnet

max_issues_repo_name: Ultramann/hnet

max_forks_repo_name: Ultramann/hnet


text:
{-# OPTIONS_GHC -fsimpl-tick-factor=150 #-}
{-# LANGUAGE BangPatterns, RecordWildCards #-}

module Network.HTTP.LoadTest.Analysis
    (
    -- * Result analysis
      Analysis(..)
    , Basic(..)
    , analyseBasic
    , analyseFull
    ) where

import Criterion.Analysis (SampleAnalysis, analyseSample)
import Network.HTTP.LoadTest.Types (Analysis(..), Basic(..), Summary(..))
import Prelude hiding (catch)
import Statistics.Quantile (weightedAvg)
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Statistics.Sample as S

analyseFull :: V.Vector Summary -> Double -> IO (Analysis SampleAnalysis)
analyseFull sumv elapsed = do
  let ci = 0.95
      resamples = 10 * 1000
  l <- analyseSample ci (G.convert . G.map summElapsed $ sumv) resamples
  return Analysis {
                 latency = l
               , latency99 = weightedAvg 99 100 . G.map summElapsed $ sumv
               , latency999 = weightedAvg 999 1000 . G.map summElapsed $ sumv
               , latValues = sumv
               , throughput = fromIntegral (G.length sumv) / elapsed
    }

analyseBasic :: V.Vector Summary -> Double -> Analysis Basic
analyseBasic sumv elapsed = Analysis {
                      latency = Basic {
                                  mean = S.mean . G.map summElapsed $ sumv
                                , stdDev = S.stdDev . G.map summElapsed $ sumv
                                }
                    , latency99 = weightedAvg 99 100 . G.map summElapsed $ sumv
                    , latency999 = weightedAvg 999 1000 . G.map summElapsed $ sumv
                    , latValues = sumv
                    , throughput = fromIntegral (G.length sumv) / elapsed
                    }

max_stars_repo_name: fhartwig/pronk

max_issues_repo_name: liqd/pronk

max_forks_repo_name: bos/pronk


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE NamedFieldPuns #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ViewPatterns #-}

module Language.Grappa.Interp.StandardHORepr where

import Data.Vector (Vector)
import qualified Data.Vector as V
import Control.Monad

import Numeric.LinearAlgebra hiding (R, Uniform, (<>), Vector)

import Language.Grappa.Distribution
import Language.Grappa.Interp
import Language.Grappa.GrappaInternals
import Language.Grappa.Frontend.DataSource

import qualified Numeric.Log as Log

import qualified Numeric.AD.Mode.Reverse as ADR
import qualified Numeric.AD.Internal.Reverse as ADR
import qualified Data.Reflection as ADR (Reifies)

----------------------------------------------------------------------
-- * The StandardHORepr Representation
----------------------------------------------------------------------

-- | The type tag for the standard higher-order representation, which is
-- parameterized by a monad and by the representation types for reals and ints
data StandardHORepr (m :: * -> *) (r :: *) (i :: *) :: *

-- | The type family for 'StandardHORepr' expressions
type family StandardHOReprF m r i a :: * where
  StandardHOReprF m r i (a -> b) =
    (StandardHOReprF m r i a -> StandardHOReprF m r i b)
  StandardHOReprF m r i (Dist a) =
    (DistVar (StandardHORepr m r i) a -> m (StandardHOReprF m r i a))
  StandardHOReprF m r i (ADT adt) =
    adt (GExpr (StandardHORepr m r i)) (ADT adt)
  StandardHOReprF m r i (Vector a) = Vector (StandardHOReprF m r i a)
  StandardHOReprF m r i Bool    = Bool
  StandardHOReprF m r i Int     = i
  StandardHOReprF m r i Prob    = Log.Log r
  StandardHOReprF m r i R       = r
  StandardHOReprF m r i RMatrix = Matrix r
  StandardHOReprF m r i a       = a

instance ValidExprRepr (StandardHORepr m r i) where
  type GExprRepr (StandardHORepr m r i) a = StandardHOReprF m r i a
  interp__'bottom = error "StandardHORepr: unexpected bottom!"
  interp__'injTuple !tup = GExpr tup
  interp__'projTuple (GExpr !tup) k = k tup
  interp__'app (GExpr !f) (GExpr x) = GExpr (f x)
  interp__'lam f = GExpr (unGExpr . f . GExpr)
  interp__'fix f = f (interp__'fix f)

instance StrongTupleRepr (StandardHORepr m r i) where
  interp__'strongProjTuple (GExpr tup) = tup

-- | Helper to match on v-expressions of atomic type in the 'StandardHORepr': if
-- the 'DistVar' is not a 'VParam', then destructure it and pass it to the
-- continuation (the 2nd argument); otherwise, return the failure continuation
-- (the 3rd argument)
matchHOReprAtomicDistVar ::
  (IsAtomic a ~ 'True, EmbedRepr (StandardHORepr m r i) a) =>
  DistVar (StandardHORepr m r i) a ->
  (GExpr (StandardHORepr m r i) a -> ret) -> ret -> ret
matchHOReprAtomicDistVar VParam _ ret = ret
matchHOReprAtomicDistVar (VData (GData a)) k _ = k $ embedRepr a
matchHOReprAtomicDistVar (VData GNoData) _ ret = ret
matchHOReprAtomicDistVar (VExpr e) k _ = k e

-- | Test if a v-expression of atomic type is a "missing value"
isMissingHOReprAtomicDistVar ::
  (IsAtomic a ~ 'True, EmbedRepr (StandardHORepr m r i) a) =>
  DistVar (StandardHORepr m r i) a -> Bool
isMissingHOReprAtomicDistVar dv =
  matchHOReprAtomicDistVar dv (\_ -> False) True

-- | Helper to match on v-expressions in the 'StandardHORepr': if the 'DistVar'
-- is not a 'VParam', then destructure it and pass it to the continuation (the
-- 2nd argument); otherwise, return the failure continuation (the 3rd argument)
matchHOReprADTDistVar ::
  TraversableADT adt =>
  DistVar (StandardHORepr m r i) (ADT adt) ->
  (adt (DistVar (StandardHORepr m r i)) (ADT adt) -> ret) -> ret -> ret
matchHOReprADTDistVar VParam _ ret = ret
matchHOReprADTDistVar (VData (GData (ADT adt))) k _ =
  k $ mapADT (VData . GData . unId) adt
matchHOReprADTDistVar (VData GNoData) _ ret = ret
matchHOReprADTDistVar (VData (GADTData adt)) k _ = k $ mapADT VData adt
matchHOReprADTDistVar (VExpr (GExpr adt)) k _ = k $ mapADT VExpr adt
matchHOReprADTDistVar (VADT adt) k _ = k adt

-- | Recursively match a v-expression with list type, returning a list of
-- 'DistVar's in the list along with a 'Bool' flag indicating whether the list
-- ends with a "missing list", i.e., a 'VParam' or @'VData' 'GNoData'@
matchHOReprListDistVar :: DistVar (StandardHORepr m r i) (GList a) ->
                          ([DistVar (StandardHORepr m r i) a], Bool)
matchHOReprListDistVar dv =
  matchHOReprADTDistVar dv
  (\adt -> case adt of
      Nil -> ([], False)
      Cons hdv tlv ->
        let (l, flag) = matchHOReprListDistVar tlv in
        (hdv:l, flag))
  ([], True)

matchHOReprVectorDistVar ::
  DistVar (StandardHORepr m r i) (Vector a) ->
  (Vector (DistVar (StandardHORepr m r i) a) -> ret) -> ret -> ret
matchHOReprVectorDistVar VParam _ ret = ret
matchHOReprVectorDistVar (VData (GData a)) k _ = k $ V.map (VData . GData) a
matchHOReprVectorDistVar (VData GNoData) _ ret = ret
matchHOReprVectorDistVar (VExpr (GExpr dvs)) k _ = k $ V.map (VExpr . GExpr) dvs

instance Monad m => ValidRepr (StandardHORepr m r i) where
  type GVExprRepr (StandardHORepr m r i) a =
    DistVar (StandardHORepr m r i) a
  type GStmtRepr (StandardHORepr m r i) a = m (StandardHOReprF m r i a)

  interp__'projTupleStmt (GExpr !tup) k = k tup

  interp__'vInjTuple !tup = GVExpr (VADT $ mapADT unGVExpr tup)
  interp__'vProjTuple (GVExpr ve) k =
    matchHOReprADTDistVar ve (k . mapADT GVExpr)
    (k $ buildTuple $ GVExpr VParam)

  interp__'vwild k = k (GVExpr VParam)
  interp__'vlift e k = k (GVExpr $ VExpr e)

  interp__'return (GExpr !x) = GStmt (return x)
  interp__'let rhs body = body rhs
  interp__'sample (GExpr !d) (GVExpr !dv) k = GStmt $ do
    !x <- d dv
    unGStmt $ k (GExpr x)

  interp__'mkDist f = GExpr (\ dv -> unGStmt $ f $ GVExpr dv)

instance TraversableADT adt =>
         Interp__ADT__Expr (StandardHORepr m r i) adt where
  interp__'injADT adt = GExpr adt
  interp__'projADT (GExpr adt) k = k adt
  interp__'projMatchADT (GExpr adt) _ matcher k_succ k_fail =
    if applyCtorMatcher matcher adt then k_succ adt else k_fail

instance (Monad m, TraversableADT adt) =>
         Interp__ADT (StandardHORepr m r i) adt where
  interp__'vInjADT adt =
    GVExpr (VADT $ mapADT unGVExpr adt)
  interp__'vProjMatchADT (GVExpr ve) ctor matcher k_succ k_fail =
    matchHOReprADTDistVar ve
    (\adt ->
      if applyCtorMatcher matcher adt then
        k_succ (mapADT GVExpr adt)
      else k_fail)
    (k_succ $ mapADT (const $ GVExpr VParam) ctor)

instance Interp__'source (StandardHORepr m Double Int) a where
  interp__'source src =
    GVExpr . VData <$> interpSource src

instance EmbedRepr (StandardHORepr m Double i) R where
  embedRepr = GExpr

instance EmbedRepr (StandardHORepr m Double i) Prob where
  embedRepr = GExpr . fromProb

instance ADR.Reifies s ADR.Tape =>
         EmbedRepr (StandardHORepr m (ADR.Reverse s Double) i) R where
  embedRepr = GExpr . ADR.auto

instance ADR.Reifies s ADR.Tape =>
         EmbedRepr (StandardHORepr m (ADR.Reverse s Double) i) Prob where
  embedRepr = GExpr . fmap ADR.auto . fromProb

instance Num i => EmbedRepr (StandardHORepr m r i) Int where
  embedRepr = GExpr . fromIntegral

instance EmbedRepr (StandardHORepr m r i) Bool where
  embedRepr = GExpr

instance MapC (EmbedRepr (StandardHORepr m r i)) ts =>
         EmbedRepr (StandardHORepr m r i) (ADT (TupleF ts)) where
  embedRepr = GExpr . helper . unADT where
    helper :: MapC (EmbedRepr (StandardHORepr m r i)) ts' =>
              TupleF ts' Id any -> TupleF ts' (GExpr (StandardHORepr m r i)) any
    helper (Tuple0) = Tuple0
    helper (Tuple1 (Id a)) = Tuple1 (embedRepr a)
    helper (Tuple2 (Id a) (Id b)) =
      Tuple2 (embedRepr a) (embedRepr b)
    helper (Tuple3 (Id a) (Id b) (Id c)) =
      Tuple3 (embedRepr a) (embedRepr b) (embedRepr c)
    helper (Tuple4 (Id a) (Id b) (Id c) (Id d)) =
      Tuple4 (embedRepr a) (embedRepr b) (embedRepr c) (embedRepr d)
    helper (TupleN (Id a) (Id b) (Id c) (Id d) (Id e) tup) =
      TupleN (embedRepr a) (embedRepr b) (embedRepr c) (embedRepr d)
      (embedRepr e) (helper tup)

instance EmbedRepr (StandardHORepr m Double i) a =>
         EmbedRepr (StandardHORepr m Double i) (Vector a) where
  embedRepr xs = GExpr $ V.map (unGExpr . helper) xs where
    helper :: EmbedRepr (StandardHORepr m Double i) a => a ->
              GExpr (StandardHORepr m Double i) a
    helper = embedRepr

----------------------------------------------------------------------
-- Boolean and comparison operations
----------------------------------------------------------------------

instance Interp__'ifThenElse (StandardHORepr m r i) where
  interp__'ifThenElse (GExpr c) t e = if c then t else e

instance Monad m => Interp__'vmatchSwitch (StandardHORepr m r Int) where
  interp__'vmatchSwitch (GExpr i) stmts = stmts !! i

instance Interp__not (StandardHORepr m r i) where
  interp__not = GExpr not

instance Interp__'amp'amp (StandardHORepr m r i) where
  interp__'amp'amp = GExpr (&&)

instance Interp__'bar'bar (StandardHORepr m r i) where
  interp__'bar'bar = GExpr (||)

instance (Eq a, Eq (StandardHOReprF m r i a)) =>
         Interp__'eq'eq (StandardHORepr m r i) a where
  interp__'eq'eq = GExpr (==)

instance (Ord a, Ord (StandardHOReprF m r i a)) =>
         Interp__'lt (StandardHORepr m r i) a where
  interp__'lt = GExpr (<)

instance (Ord a, Ord (StandardHOReprF m r i a)) =>
         Interp__'gt (StandardHORepr m r i) a where
  interp__'gt = GExpr (>)

instance (Ord a, Ord (StandardHOReprF m r i a)) =>
         Interp__'lt'eq (StandardHORepr m r i) a where
  interp__'lt'eq = GExpr (<=)

instance (Ord a, Ord (StandardHOReprF m r i a)) =>
         Interp__'gt'eq (StandardHORepr m r i) a where
  interp__'gt'eq = GExpr (>=)

instance (Ord a, Ord (StandardHOReprF m r i a)) =>
         Interp__min (StandardHORepr m r i) a where
  interp__min = GExpr min

instance (Ord a, Ord (StandardHOReprF m r i a)) =>
         Interp__max (StandardHORepr m r i) a where
  interp__max = GExpr max

----------------------------------------------------------------------
-- Numeric Operations
----------------------------------------------------------------------

instance (Num a, Num (StandardHOReprF m r i a)) =>
         Interp__'plus (StandardHORepr m r i) a where
  interp__'plus = GExpr (+)

instance (Num a, Num (StandardHOReprF m r i a)) =>
         Interp__'minus (StandardHORepr m r i) a where
  interp__'minus = GExpr (-)

instance (Num a, Num (StandardHOReprF m r i a)) =>
         Interp__'times (StandardHORepr m r i) a where
  interp__'times = GExpr (*)

instance (Num a, Num (StandardHOReprF m r i a)) =>
         Interp__negate (StandardHORepr m r i) a where
  interp__negate = GExpr negate

instance (Num a, Num (StandardHOReprF m r i a)) =>
         Interp__abs (StandardHORepr m r i) a where
  interp__abs = GExpr abs

instance (Num a, Num (StandardHOReprF m r i a)) =>
         Interp__signum (StandardHORepr m r i) a where
  interp__signum = GExpr signum

instance (Num a, Num (StandardHOReprF m r i a)) =>
         Interp__fromInteger (StandardHORepr m r i) a where
  interp__fromInteger = GExpr fromInteger

instance (Num a, Num (StandardHOReprF m r i a)) =>
         Interp__'integer (StandardHORepr m r i) a where
  interp__'integer n = GExpr (fromInteger n)

instance (Interp__'integer (StandardHORepr m r i) a,
          Eq (StandardHOReprF m r i a))
         => Interp__'eqInteger (StandardHORepr m r i) a where
  interp__'eqInteger (GExpr x) (GExpr y) = GExpr (x == y)

instance (Fractional a, Fractional (StandardHOReprF m r i a)) =>
         Interp__'div (StandardHORepr m r i) a where
  interp__'div = GExpr (/)

instance (Fractional a, Fractional (StandardHOReprF m r i a)) =>
         Interp__recip (StandardHORepr m r i) a where
  interp__recip = GExpr recip

instance (Fractional a, Fractional (StandardHOReprF m r i a)) =>
         Interp__fromRational (StandardHORepr m r i) a where
  interp__fromRational = GExpr fromRational

instance (Fractional a, Fractional (StandardHOReprF m r i a)) =>
         Interp__'rational (StandardHORepr m r i) a where
  interp__'rational n = GExpr (fromRational n)

instance (Interp__'rational (StandardHORepr m r i) a,
          Eq (StandardHOReprF m r i a))
         => Interp__'eqRational (StandardHORepr m r i) a where
  interp__'eqRational (GExpr x) (GExpr y) = GExpr (x == y)

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__pi (StandardHORepr m r i) a where
  interp__pi = GExpr pi

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__exp (StandardHORepr m r i) a where
  interp__exp = GExpr exp

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__log (StandardHORepr m r i) a where
  interp__log = GExpr log

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__sqrt (StandardHORepr m r i) a where
  interp__sqrt = GExpr sqrt

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__'times'times (StandardHORepr m r i) a where
  interp__'times'times = GExpr (**)

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__logBase (StandardHORepr m r i) a where
  interp__logBase = GExpr logBase

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__sin (StandardHORepr m r i) a where
  interp__sin = GExpr sin

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__cos (StandardHORepr m r i) a where
  interp__cos = GExpr cos

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__tan (StandardHORepr m r i) a where
  interp__tan = GExpr tan

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__asin (StandardHORepr m r i) a where
  interp__asin = GExpr asin

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__acos (StandardHORepr m r i) a where
  interp__acos = GExpr acos

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__atan (StandardHORepr m r i) a where
  interp__atan = GExpr atan

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__sinh (StandardHORepr m r i) a where
  interp__sinh = GExpr sinh

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__cosh (StandardHORepr m r i) a where
  interp__cosh = GExpr cosh

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__tanh (StandardHORepr m r i) a where
  interp__tanh = GExpr tanh

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__asinh (StandardHORepr m r i) a where
  interp__asinh = GExpr asinh

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__acosh (StandardHORepr m r i) a where
  interp__acosh = GExpr acosh

instance (Floating a, Floating (StandardHOReprF m r i a)) =>
         Interp__atanh (StandardHORepr m r i) a where
  interp__atanh = GExpr atanh

----------------------------------------------------------------------
-- Probability expression operations
----------------------------------------------------------------------

instance (Ord r, Floating r) => Interp__realToProb (StandardHORepr m r i) where
  interp__realToProb = GExpr (Log.Exp . log . toNonNeg) where
    toNonNeg :: (Ord r, Floating r) => r -> r
    toNonNeg r = if r < 0 then 0 else r

instance Interp__logRealToProb (StandardHORepr m r i) where
  interp__logRealToProb = GExpr Log.Exp

instance Floating r => Interp__probToReal (StandardHORepr m r i) where
  interp__probToReal = GExpr (exp . Log.ln)

instance Interp__probToLogReal (StandardHORepr m r i) where
  interp__probToLogReal = GExpr Log.ln

instance HasGamma r => Interp__gammaProb (StandardHORepr m r i) where
  interp__gammaProb = GExpr (Log.Exp . logGamma)

instance HasGamma r => Interp__digamma (StandardHORepr m r i) where
  interp__digamma = GExpr digamma

----------------------------------------------------------------------
-- Misc operations
----------------------------------------------------------------------

instance (Show a, Show (StandardHOReprF m r i a), i ~ Int) =>
         Interp__gtrace (StandardHORepr m r i) a b where
  interp__gtrace = GExpr gtrace

instance i ~ Int => Interp__gerror (StandardHORepr m r i) a where
  interp__gerror = GExpr gerror

----------------------------------------------------------------------
-- * Distributions Supported by All 'StandardHORepr' Representations
----------------------------------------------------------------------

instance Monad m => Interp__ctorDist__ListF (StandardHORepr m r i) where
  interp__ctorDist__Nil = GExpr $ \ mkNil dv ->
    do _ <-
         matchHOReprADTDistVar dv
         (\adt -> case adt of
             Nil -> mkNil (VADT Tuple0)
             Cons _ _ -> error "Unexpected Cons")
         (mkNil VParam)
       return Nil

  interp__ctorDist__Cons = GExpr $ \ mkCons dv ->
    do (Tuple2 hd tl) <-
         matchHOReprADTDistVar dv
         (\adt -> case adt of
             Nil -> error "Unexpected Nil"
             Cons hdv tlv -> mkCons (VADT (Tuple2 hdv tlv)))
         (mkCons (VADT (Tuple2 VParam VParam)))
       return (Cons hd tl)

-- If a repr can do a categorical, it can do an ADT distribution on lists
instance (Monad m, Num i, Eq i, Show i,
          Interp__categorical (StandardHORepr m r i)) =>
         Interp__adtDist__ListF (StandardHORepr m r i) where
  interp__adtDist__ListF =
    GExpr $ \ probNil mkNil probCons mkCons dvList ->
    let
      -- Build a categorical distribution for the constructor, where 0 -> Nil
      -- and 1 -> Cons
      ctor_dist :: DistVar (StandardHORepr m r i) Int ->
                   m (StandardHOReprF m r i Int)
      ctor_dist =
        unGExpr (interp__categorical
                 :: GExpr (StandardHORepr m r i) (GList Prob -> Dist Int)) $
        fromHaskellListF GExpr [GExpr probNil, GExpr probCons]
      -- Helper wrapper around mkNil
      mkNilH = mkNil (VADT Tuple0)
      -- Helper wrapper around mkCons, that takes vars for the head and tail
      mkConsH hdv tlv =
        do Tuple2 hd tl <- mkCons (VADT (Tuple2 hdv tlv))
           return (Cons hd tl) in

    matchHOReprADTDistVar dvList
    (\adt -> case adt of
        Nil ->
          void (ctor_dist $ VData $ GData 0) >> void mkNilH >> return Nil
        Cons hdv tlv ->
          void (ctor_dist $ VData $ GData 1) >> mkConsH hdv tlv)
    (do ctor_choice <- ctor_dist VParam
        case ctor_choice of
          0 -> void mkNilH >> return Nil
          1 -> mkConsH VParam VParam
          _ -> error ("ListF: Invalid constructor choice: "
                      ++ show ctor_choice))

--
-- * Vector Operations
--

instance (Monad m, i ~ Int, EmbedRepr (StandardHORepr m r Int) a) =>
         Interp__vec_iid (StandardHORepr m r i) a where
  interp__vec_iid = GExpr $ \len d dv ->
    matchHOReprVectorDistVar dv
    (\dvs ->
      if V.length dvs == len then V.mapM d dvs else
        error "vec_iid: incorrect vector length on input")
    (V.replicateM len (d VParam))

{-
-- | Build a distribution on 'Vector's from one on lists
instance (Monad m, EmbedRepr (StandardHORepr m r Int) a) =>
         Interp__vec_dist (StandardHORepr m r Int) a where
  interp__vec_dist = GExpr $ \d dv ->
    matchHOReprVectorDistVar dv
    (\dvs ->
      V.fromList <$> toHaskellListF unGExpr <$> unGExpr <$>
      d _)
    _
-}

instance Monad m => Interp__vec_nil_dist (StandardHORepr m r i) a where
  interp__vec_nil_dist = GExpr $ \dv ->
    matchHOReprVectorDistVar dv
    (\dvs ->
      if V.length dvs == 0 then return V.empty else
        error "vec_nil_dist: non-empty vector!")
    (return V.empty)

instance Monad m => Interp__vec_cons_dist (StandardHORepr m r i) a where
  interp__vec_cons_dist = GExpr $ \d dv ->
    do (Tuple2 (GExpr hd) (GExpr tl)) <-
         case dv of
           VParam -> d VParam
           VData (GData xs) ->
             if V.length xs == 0 then
               error "vec_cons_dist: empty vector!"
             else
               d $ VData $ GData $ ADT $ Tuple2 (Id $ V.head xs) (Id $ V.tail xs)
           VData GNoData -> d (VData GNoData)
           VExpr (GExpr es) ->
             if V.length es == 0 then
               error "vec_cons_dist: empty vector!"
             else
               d $ VExpr $ GExpr $ Tuple2 (GExpr $ V.head es) (GExpr $ V.tail es)
       return (V.cons hd tl)

instance Interp__vec_head (StandardHORepr m r i) a where
  interp__vec_head = GExpr V.head

instance Interp__vec_tail (StandardHORepr m r i) a where
  interp__vec_tail = GExpr V.tail

instance i ~ Int => Interp__vec_length (StandardHORepr m r i) a where
  interp__vec_length = GExpr V.length

--
-- * Matrix Operations
--

-- FIXME: matrix operations!

max_stars_repo_name: kquick/grappa

max_issues_repo_name: kquick/grappa

max_forks_repo_name: kquick/grappa


text:
module BuckinghamPi
    (
    FundamentalUnit(..),
    DimensionalVariable(..),
    DimensionlessVariable(..),
    Exponent(..),
    generatePiGroups,
    numFundamentalUnits
    ) where

import           Control.Applicative
import           Data.List                     (nub, (\\))
import           Numeric.LinearAlgebra.Data
import           Numeric.LinearAlgebra.HMatrix
import           Text.Printf                   (printf)

-- Define sum type for fundamental units (M,L,T,Theta).
data FundamentalUnit = Mass | Length | Time | Temperature deriving (Eq,Show)

-- Define type synonym for Exponent.
type Exponent = Double

-- A variable is defined by the name, along with units, which is a list of tuples which contains
-- the FundamentalUnit along with an exponent.
data DimensionalVariable = DimensionalVariable {name :: String, units :: [(FundamentalUnit, Exponent)]} deriving (Eq,Show)

-- Define data type for a dimensionless variable.
data DimensionlessVariable = DimensionlessVariable {dimensionalVariables :: [DimensionalVariable], exponents :: [Exponent]} deriving (Eq)

-- Define function which normalizes exponents so that the smallest value is +/- 1.
normalizeExponents :: [Exponent] -> [Exponent]
normalizeExponents exps = map (/factor) exps
  where
    minExp = minimum exps
    factor = minExp * signum minExp

-- Define show instance for a dimensionless variable.
instance Show DimensionlessVariable where
  show (DimensionlessVariable dimVars exps) = concatMap (uncurry (printf "%s^(%.4f) ")) $ filter (\t -> abs (snd t) > 1.0e-6) $ zip (map name dimVars) (normalizeExponents exps)

getFundamentalUnits :: DimensionalVariable -> [FundamentalUnit]
getFundamentalUnits dimVar = map fst $ units dimVar

getUniqueFundamentalUnits :: [DimensionalVariable] -> [FundamentalUnit]
getUniqueFundamentalUnits = nub . concatMap getFundamentalUnits

numFundamentalUnits :: [DimensionalVariable] -> Int
numFundamentalUnits = length . getUniqueFundamentalUnits

-- Define function that computes the number of dimensionless groups given a list of dimensional variables.
computeNumberOfGroups :: [DimensionalVariable] -> Int
computeNumberOfGroups dimensionalVars = length dimensionalVars - numFundamentalUnits dimensionalVars

-- Define function that returns nonrepeating variables.
getNonrepeatingVars :: [DimensionalVariable] -> [Int] -> [DimensionalVariable]
getNonrepeatingVars dimVars repeatedIndices = [dimVars !! i | i <- [0..(length dimVars - 1)], i `notElem` repeatedIndices]

-- Define a function that returns the exponent on a given unit in the given dimensional variable.
getExponent :: FundamentalUnit -> DimensionalVariable -> Exponent
getExponent unit dimVar
  | unit `elem` map fst (units dimVar) = snd $ head $ filter (\u -> fst u == unit) $ units dimVar
  | otherwise = 0.0

-- Define a function that returns the equation that must be solved.
getEquation ::  [DimensionalVariable] -> FundamentalUnit -> [Exponent]
getEquation varCombination funUnit = map (getExponent funUnit) varCombination

getVariableCombinations :: [DimensionalVariable] -> [DimensionalVariable] -> [[DimensionalVariable]]
getVariableCombinations repeatingVars = map (\x -> repeatingVars ++ [x])

buildMatrix :: [DimensionalVariable] -> [FundamentalUnit] -> Matrix Double
buildMatrix varCombination funUnits = fromLists $ map (getEquation varCombination) funUnits

findDimensionlessGroup :: [FundamentalUnit] -> [DimensionalVariable] -> DimensionlessVariable
findDimensionlessGroup funUnits dimVars = DimensionlessVariable dimVars (concat . toLists $ nullspace $ buildMatrix dimVars funUnits)

-- Define function which takes a list of dimensional variables along with a list of integers
-- that identifies which are the repeated variables, and returns a list of dimensionless variables.
generatePiGroups :: [DimensionalVariable] -> [Int] -> [DimensionlessVariable]
generatePiGroups dimVars repeatingVarInds = map (findDimensionlessGroup funUnits) (getVariableCombinations repeatingVars nonrepeatingVars)
  where
    repeatingVars = [dimVars !! i | i <- repeatingVarInds]
    nonrepeatingVars = dimVars \\ repeatingVars
    funUnits = getUniqueFundamentalUnits dimVars


text:
module Main where

import Vector
import Tensor
import Data.Complex
import System.Random.Stateful
import EitherTrans
import Control.Monad.Trans.State

func i n c = do
        s <- newIOGenM $ mkStdGen 42
        let sT = measure i s
        result <- runEither $ (EitherT . return) x1 >>= execStateT sT
        let result2 = x1 >>= isEntangled [c]
        let result3 = x2 >>= isEntangled [c]
        --either putStrLn (putStrLn . illegalPeek) result3
        either putStrLn print result2
        either putStrLn print result3
    where
        x = initNumQubit0 n :: Either String (Vector (Complex Double))
        x0 = x >>= applyGate hadamard 0
        x1 = x0 >>= applyControl pauliX 0 1
        x2 = x >>= applyGateAll hadamard

main = func 1 2 0


text:
module Main where

import Neural.Matrix
import Numeric.LinearAlgebra
import Test.Hspec

{- Simple Net
   1
      [3, 2, -2] -> 1  |-1.5| = -0.5 -> 0                 bias
   2                    bias              [-2, 0.5] -> 7 | -3 | -> 4
      [-1, 2, 3] -> 12 | +2 | =  14 -> 14
   3
-}

zeroVec = vector [0.0,0.0]
zeroNet = [(relu, zeroVec, matrix 2 [0,0,0,0]), (relu, vector [0], matrix 1 [0,0])]
vTest = vector [1, 2, 3]
mTest = matrix 2 [3.0, -1.0, 2.0, 2.0, -2.0, 3.0]
simpleNet = [(relu, vector [-1.5, 2.0], mTest), (relu, vector [-3.0], matrix 1 [-2.0, 0.5])]

main :: IO ()
main = hspec $ do
  describe "Check feed layer" $ do
    it "passes trivial case" $ do
      (feedLayer (zeroVec, zeroVec) (head zeroNet)) `shouldBe` (zeroVec, zeroVec)
    it "passes simple case" $ do
      (feedLayer (vTest, vTest) (head simpleNet)) `shouldBe` (vector [-0.5,14.0], vector [0.0,14.0])

  describe "Check forward pass" $ do
    it "passes trivial case" $ do
      let zp = [(zeroVec, zeroVec), (zeroVec, zeroVec),
                (vector [0.0], vector [0.0])]
      (forwardPass zeroVec zeroNet) `shouldBe` zp
    it "passes simple case" $ do
      let fp = [(vector [1.0,2.0,3.0], vector [1.0,2.0,3.0]),
                (vector [-0.5,14.0], vector [0.0,14.0]),
                (vector [4.0], vector [4.0])]
      (forwardPass vTest simpleNet) `shouldBe` fp

  describe "Check chain" $ do
    it "passes trivial case" $ do
      let (_, _, zm) = last zeroNet
      let ze = vector [0.0]
      (chain (zeroVec, zm) ze) `shouldBe` zeroVec
    it "passes simple case" $ do
      let tl = vector [-0.5,14]
      let (_, _, lw) = last simpleNet
      let e = vector [0.5]
      let chained = vector [0,0.25]
      (chain (tl, lw) e) `shouldBe` chained

  describe "Check back prop" $ do
    it "passes simple case" $ do
      let tls = [zeroVec, vector [0]]
      let lws = map (\(_, _, x) -> x) zeroNet
      (backProp (vector [0]) tls lws) `shouldBe` tls


text:
-- File created: 2009-07-11 21:47:19

-- ScmValue depends on Haschoo
-- Haschoo  depends on Context
-- Context  depends on ScmValue
--
-- Thus they all have to be in the same module, to avoid circular dependencies.
{-# LANGUAGE BangPatterns, GeneralizedNewtypeDeriving #-}

module Haschoo.Types
   ( Haschoo, runHaschoo, withHaschoo
   , ScmValue(..), MacroCall(..)
   , isTrue, isAggregate
   , toScmString, toScmMString
   , toScmVector
   , pairToList, listToPair
   , Context, mkContext, addToContext, contextLookup, contextSize
   , scmShow, scmShowWith
   ) where

import Control.Applicative        ((<$>))
import Control.Arrow              (second)
import Control.Monad.Error        (ErrorT, MonadError, runErrorT)
import Control.Monad.State.Strict (StateT, MonadState, evalStateT)
import Control.Monad.Trans        (MonadIO, liftIO)
import Data.Array.IO              (IOArray, IOUArray, newListArray, getElems)
import Data.Array.Unboxed         (Array, UArray, elems, listArray)
import Data.Complex               (Complex((:+)))
import Data.IORef                 (IORef, readIORef, newIORef)
import Data.List                  (intercalate)
import Data.Ratio                 (numerator, denominator)
import System.IO                  (Handle)
import System.IO.Unsafe           (unsafeInterleaveIO)
import Text.Show.Functions        ()

import qualified Data.ListTrie.Patricia.Map.Enum as TM
import Data.ListTrie.Patricia.Map.Enum (TrieMap)

import Haschoo.Utils (ErrOr, lazyMapM)

-- Our main monad
newtype Haschoo a = Haschoo (StateT HaschState (ErrorT String IO) a)
   deriving (Functor, Monad, MonadIO, MonadState HaschState, MonadError String)

runHaschoo :: HaschState -> Haschoo a -> IO (ErrOr a)
runHaschoo state (Haschoo h) = runErrorT $ evalStateT h state

-- Like 'local' of the Reader monad
withHaschoo :: HaschState -> Haschoo a -> Haschoo a
withHaschoo state h = either fail return =<< liftIO (runHaschoo state h)

type HaschState = [IORef Context]

data ScmValue = ScmPrim  String !([ScmValue] -> Haschoo   ScmValue)
              | ScmFunc  String !([ScmValue] -> IO (ErrOr ScmValue))

              | ScmMacro String
                   !(MacroCall ->
                        Haschoo (ScmValue, TrieMap Char [IORef Context]))

              | ScmBool       !Bool
              | ScmChar       !Char
              | ScmInt        !Integer
              | ScmRat        !Rational
              | ScmReal       !Double
              | ScmComplex    !(Complex Double)
              | ScmIdentifier !String
              | ScmPair       !(IORef ScmValue) !(IORef ScmValue)
              | ScmString     !(  UArray Int Char)
              | ScmMString    !(IOUArray Int Char)
              | ScmVector     !(   Array Int ScmValue)
              | ScmMVector    !( IOArray Int ScmValue)

              | ScmSyntax { patterns :: ![(ScmValue, ScmValue, [String])]
                          , literals :: ![(String, Maybe ScmValue)] }

              -- These two are only for literal lists, appearing only directly
              -- from the parser
              | ScmList       ![ScmValue]
              | ScmDottedList ![ScmValue] !ScmValue -- The list is never empty

              -- The "unspecified value" returned by IO procedures and such
              | ScmVoid

              | ScmInput  !Handle
              | ScmOutput !Handle
              | ScmEOF

              -- The "environment-specifier" used by eval-related procedures
              | ScmContext ![Context]

-- The form of the macro usage:
--
-- (macro foo bar), (macro foo . bar), or #(macro foo bar)
--
-- The macro keyword is not included regardless of constructor
data MacroCall = MCList   ![ScmValue]
               | MCDotted ![ScmValue] !ScmValue

isTrue, isAggregate :: ScmValue -> Bool
isTrue (ScmBool False) = False
isTrue _               = True

isAggregate (ScmList       _)   = True
isAggregate (ScmDottedList _ _) = True
isAggregate (ScmVector     _)   = True
isAggregate (ScmMVector    _)   = True
isAggregate (ScmPair       _ _) = True
isAggregate _                   = False

-- From ScmPair to Left . ScmDottedList or Right . ScmList
pairToList :: ScmValue -> IO (Either ScmValue ScmValue)
pairToList = go []
 where
   go xs (ScmList       ys)   =
      return . Right $ ScmList       (reverse xs ++ ys)
   go xs (ScmDottedList ys z) =
      return . Left  $ ScmDottedList (reverse xs ++ ys) z

   go xs (ScmPair a b) = do
      a' <- readIORef a
      b' <- readIORef b
      go (a':xs) b'

   go xs v = return . Left $ ScmDottedList (reverse xs) v

-- ScmList -> ScmPair, returns the tail pointer if it wasn't empty
listToPair :: [ScmValue] -> IO (ScmValue, Maybe (IORef ScmValue))
listToPair []     = return (ScmList [], Nothing)
listToPair (x:xs) = second Just <$> go x xs
 where
   go :: ScmValue -> [ScmValue] -> IO (ScmValue, IORef ScmValue)
   go v [] = do
      a <- newIORef v
      b <- newIORef (ScmList [])
      return (ScmPair a b, b)

   go v (y:ys) = do
      a         <- newIORef v
      (zs, fin) <- go y ys
      b         <- newIORef zs
      return (ScmPair a b, fin)

toScmString :: String -> ScmValue
toScmString s = ScmString $ listArray (0, length s - 1) s

toScmMString :: String -> IO ScmValue
toScmMString s = ScmMString <$> newListArray (0, length s - 1) s

toScmVector :: [ScmValue] -> ScmValue
toScmVector v = ScmVector $ listArray (0, length v - 1) v

newtype Context = Context (TrieMap Char ScmValue)

mkContext :: [(String, ScmValue)] -> Context
mkContext = Context . TM.fromList

addToContext :: String -> ScmValue -> Context -> Context
addToContext name val (Context c) = Context (TM.insert' name val c)

contextLookup :: String -> Context -> Maybe ScmValue
contextLookup s (Context c) = TM.lookup s c

contextSize :: Context -> Int
contextSize (Context c) = TM.size' c

scmShow :: ScmValue -> IO String
scmShow ScmVoid           = return "" -- Has no representation
scmShow (ScmBool b)       = return$ '#' : if b then "t" else "f"
scmShow (ScmPrim s _)     = return s
scmShow (ScmFunc s _)     = return s
scmShow (ScmMacro s _)    = return s
scmShow (ScmIdentifier s) = return s
scmShow (ScmInt n)        = return$ show n
scmShow (ScmReal n)
   | isNaN      n = return "+nan.#"
   | isInfinite n = return$ (if n < 0 then '-' else '+') : "inf.#"
   | otherwise    = return$ show n
scmShow (ScmRat n) =
   if denominator n == 1
      then return$ show (numerator n)
      else return$ concat [show (numerator n), "/", show (denominator n)]

scmShow (ScmComplex (a :+ b)) = do
   bs <- scmShow (ScmReal b)
   as <- scmShow (ScmReal a)
   return$ concat [ as
                  , if head bs `elem` "-+" then [] else "+"
                  , bs
                  , "i" ]

scmShow (ScmChar c) | c == ' '  = return  "#\\space"
                    | c == '\n' = return  "#\\newline"
                    | otherwise = return$ "#\\" ++ [c]

scmShow (ScmSyntax  _ _) = return "<syntax rules>"
scmShow (ScmContext _)   = return "<environment specifier>"
scmShow ScmEOF           = return "<end of file object>"
scmShow (ScmInput  _)    = return "<input port>"
scmShow (ScmOutput _)    = return "<output port>"

scmShow (ScmString  s) = return$ showScmString (elems s)
scmShow (ScmMString s) = showScmString <$> getElems s

scmShow x | isAggregate x = scmShowWith scmShow x
          | otherwise     = error "scmShow :: the impossible happened"

scmShowWith :: (ScmValue -> IO String) -> ScmValue -> IO String
scmShowWith f (ScmList xs)        = showScmList f xs
scmShowWith f (ScmDottedList a b) = do
   as <- showScmList f a
   bs <- f b
   return$ concat [init as, " . ", bs, ")"]

scmShowWith f (ScmVector  v) = ('#':) <$>  showScmList f       (elems v)
scmShowWith f (ScmMVector v) = ('#':) <$> (showScmList f =<< getElems v)

scmShowWith f (ScmPair car cdr) = do
   a  <- readIORef car
   as <- f a
   ('(':) . (as ++) <$> unsafeInterleaveIO (readIORef cdr >>= go)
 where
   go   (ScmList [])        = return ")"
   go x@(ScmList _)         = (' ':) . tail <$> f x
   go x@(ScmDottedList _ _) = (' ':) . tail <$> f x

   go (ScmPair x y) = do
      a <- readIORef x
      as <- f a
      (' ':) . (as ++) <$> unsafeInterleaveIO (readIORef y >>= go)

   go x = (" . "++) . (++")") <$> f x

scmShowWith _ _ = error "scmShowWith :: the impossible happened"

showScmString :: String -> String
showScmString s = '"' : foldr ((.) . f) id s "\""
 where
   f c | c == '\\' = showString "\\\\"
       | c == '\"' = showString "\\\""
       | otherwise = showChar c

showScmList :: (a -> IO String) -> [a] -> IO String
showScmList f xs =
   ("("++) . (++")") . intercalate " " <$> lazyMapM f xs

max_stars_repo_name: Deewiant/haschoo

max_issues_repo_name: Deewiant/haschoo

max_forks_repo_name: Deewiant/haschoo


text:
{-# language LambdaCase #-}
{-# language FlexibleInstances #-}
{-# language OverloadedStrings #-}
{-# language GeneralizedNewtypeDeriving #-}
{-# language TypeFamilies #-}
{-# language DerivingVia #-}
{-# language DeriveTraversable #-}
{-# language UndecidableInstances #-}
{-# language FlexibleContexts #-}

module Geometry.Decorations
  ( WithOptions (..)
  , Options, Option(..)
  , mkOptions, getOptions
  , Decorated(..), fromDecorated
  , Decorator(..)
  , (#:), (#::)
  , visible, invisible
  , stroke, white, fill
  , thickness, thin
  , dashed, dotted
  , arcs
  , label, loffs, lpos, lparam
  )
where

import Data.Monoid
import Data.Complex
import Data.Maybe
import Data.String( IsString(..) )
import Control.Monad
import Data.Functor.Identity

import Geometry.Base
import Geometry.Point
import Geometry.Line
import Geometry.Circle
import Geometry.Angle
import Geometry.Polygon
import Geometry.Plot

-- | Possible SVG options for a figure.
data Option = Invisible Bool
            | Stroke String
            | Fill String
            | Thickness String
            | Dashing String
            | MultiStroke Int
            | LabelText String
            | LabelCorner (Int, Int)
            | LabelPosition Cmp
            | LabelOffset Cmp
            | LabelAngle Direction
            | SegmentMark Int
            deriving (Show)

-- | Monoidal options list wrapper. Concatenates dually to a usual list.
newtype Options = Options (Dual [Option])
  deriving (Semigroup, Monoid, Show)

-- | Puts a list of options to a wrapper.g
mkOptions = Options . Dual
-- | Extracts an option list from the `Options` wrapper.
getOptions (Options (Dual os)) = os

-- | The class of objects which could have decorations.
class WithOptions a where
  -- | Get the decoration data
  options :: a -> Options
  options = defaultOptions

  -- | Set the decoration data
  setOptions :: Options -> a -> a
  setOptions _ = id

  -- | Default labeling settings
  defaultOptions :: a -> Options
  defaultOptions _ = mempty

instance
  {-# OVERLAPPABLE #-}
  (Functor t, Foldable t, WithOptions a) => WithOptions (t a) where
  options = foldMap options
  setOptions o a = setOptions o <$> a

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

-- | The transparent decoration wrapper for geometric objects.
-- Inherits all properties of embedded object.
newtype Decorated a = Decorated (Options, a)
  deriving (Functor, Foldable, Traversable)

-- | The selector for the embedded object.
fromDecorated (Decorated (_, x)) = x

instance Applicative Decorated where
  pure p = Decorated (mempty, p)
  (<*>) = ap

instance Monad Decorated where
  Decorated (d, x) >>= f =
    let Decorated (d', y) = f x
    in Decorated (d <> d', y)

instance {-# OVERLAPPING #-} WithOptions (Decorated a) where
  options (Decorated (o, _)) = o
  setOptions o' f = Decorated (o', id) <*> f

instance Show a => Show (Decorated a) where
  show f = l <> show (fromDecorated f)
    where l = maybe mempty (<> ":") lab
          lab = find optLabelText f
          find p d = extractOption p $ defaultOptions d <> options d
          extractOption p = getFirst . foldMap (First . p) . getOptions
          optLabelText = \case {LabelText x -> Just x; _ -> Nothing }

instance Eq a => Eq (Decorated a) where
  d1 == d2 = fromDecorated d1 == fromDecorated d2

instance Metric a => Metric (Decorated a) where
  dist a b = dist (fromDecorated a) (fromDecorated b)
  dist2 a b = dist2 (fromDecorated a) (fromDecorated b)

instance Affine a => Affine (Decorated a) where
  cmp = cmp . fromDecorated
  asCmp = pure . asCmp

instance {-# OVERLAPPING #-}  Trans a => Trans (Decorated a) where
  transform t = fmap (transform t)

instance Manifold m => Manifold (Decorated m) where
  type Domain (Decorated m) = Domain m
  param = param . fromDecorated
  project = project . fromDecorated
  paramMaybe = paramMaybe . fromDecorated
  projectMaybe = projectMaybe . fromDecorated
  isContaining = isContaining . fromDecorated
  unit = unit . fromDecorated

instance Curve c => Curve (Decorated c) where
  tangent = tangent . fromDecorated
  normal = normal . fromDecorated

instance ClosedCurve c => ClosedCurve (Decorated c) where
  isEnclosing = isEnclosing . fromDecorated
  location = location . fromDecorated

instance {-# OVERLAPPING #-} Figure a => Figure (Decorated a) where
  isTrivial = isTrivial . fromDecorated
  refPoint = refPoint . fromDecorated
  box = box . fromDecorated

instance APoint p => APoint (Decorated p) where
  toPoint = toPoint . fromDecorated
  asPoint = pure . asPoint

instance Linear l => Linear (Decorated l) where
  refPoints = refPoints . fromDecorated

instance Circular a => Circular (Decorated a) where
  toCircle = toCircle . fromDecorated
  asCircle = pure . asCircle

instance PiecewiseLinear a => PiecewiseLinear (Decorated a) where
  vertices = vertices . fromDecorated
  asPolyline = asPolyline . fromDecorated

instance Polygonal p => Polygonal (Decorated p) where

instance Angular a => Angular (Decorated a) where
  asAngle = pure . asAngle
  toAngle = toAngle . fromDecorated
  setValue v = fmap (setValue v)

instance APlot p => APlot (Decorated p) where
  rmap = fmap . rmap

------------------------------------------------------------
-- | A wrapped decoration function with monoidal properties,
-- corresponding to decoration options.
newtype Decorator a = Decorator (a -> Decorated a)

mkDecorator opt val = Decorator $
  \d -> setOptions (mkOptions [ opt val ]) (pure d)

instance Semigroup (Decorator a) where
  Decorator a <> Decorator b = Decorator (a >=> b)

instance Monoid (Decorator a) where
  mempty = Decorator $ \a -> Decorated (mempty, a)

instance WithOptions a => IsString (Decorator a) where
  fromString = label

infixl 5 #:
{- | The infix operator for decorator application.

>>> aPoint # at (4, 5) #: label "A"
A:<Point (4.0, 5.0)>

>>> segment (4,5) (6,9) #: "s" <> dotted <> white
s:<Segment (4.0,5.0) (6.0,9.0)>
-}
(#:) :: WithOptions a => a -> Decorator a -> Decorated a
a #: (Decorator d) = d a

(#::) :: Decorated a -> Decorator (Decorated a) -> Decorated a
a #:: (Decorator d) = join (d a)

-- | The stroke color decorator.
stroke :: WithOptions a => String -> Decorator a
stroke = mkDecorator Stroke

-- | The decorator for white lines.
white :: WithOptions a => Decorator a
white = stroke "white"

-- | The fill color decorator.
fill :: WithOptions a => String -> Decorator a
fill = mkDecorator Fill

-- | The stroke-thickness decorator.
thickness :: WithOptions a => String -> Decorator a
thickness = mkDecorator Thickness

-- | The decorator for thin lines.
thin :: WithOptions a => Decorator a
thin = thickness "1"

-- | The decorator for dashed lines.
dashed :: WithOptions a => Decorator a
dashed = mkDecorator Dashing "5,5"

-- | The decorator for dotted lines.
dotted :: WithOptions a => Decorator a
dotted = mkDecorator Dashing "2,3"

-- | The decorator for dotted lines.
arcs :: Int -> Decorator Angle
arcs = mkDecorator MultiStroke

-- | The decorator which makes invisible object visible.
visible :: WithOptions a => Decorator a
visible = mkDecorator Invisible False

-- | The decorator which makes visible object invisible.
invisible :: WithOptions a => Decorator a
invisible = mkDecorator Invisible True

-- | The decorator for labeling objects. Could be used as overloaded string.
--
-- >>> aPoint # at (4, 5) #: "A"
-- A:<Point (4, 5)>
--
label :: WithOptions a => String -> Decorator a
label = mkDecorator LabelText

-- | The decorator for label offset.
loffs :: (Affine p,  WithOptions a) => p -> Decorator a
loffs = mkDecorator LabelOffset . conjugate . cmp

-- | The decorator for label position.
lpos :: (Affine p, WithOptions a) => p -> Decorator a
lpos = mkDecorator LabelPosition . cmp

-- | The decorator for setting label on a curve at a given parameter value.
lparam :: (WithOptions m, Manifold m) => Double -> Decorator m
lparam x = Decorator $ \f -> f #: lpos (f @-> x)

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

max_stars_repo_name: MethaHardworker/geometry

max_issues_repo_name: MethaHardworker/geometry

max_forks_repo_name: MethaHardworker/geometry


text:
------------------------------------------------
-- |
-- Module    :  Numeric.Matrix.Integral
-- Copyright :  (c) Jun Yoshida 2019
-- License   :  BSD3
--
-- Operations on integer matrices
--
------------------------------------------------

module Numeric.Matrix.Integral (
  -- Hermite normal form
  hermiteNF,
  -- Smith normal form
  smithNF,
  smithRep,
  ) where

import Control.Monad
import Control.Monad.ST (ST, runST)
import Control.Monad.Loops (whileM_)
import Data.STRef

import qualified Numeric.LinearAlgebra as LA

--import Numeric.Matrix.Integral.HNFLLL
--import Numeric.Matrix.Integral.SmithNF
import Numeric.Matrix.Integral.NormalForms

-- 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 LA.! d LA.! d) /= 0
  whileM_ (p <$> readSTRef dRef) $ modifySTRef' dRef (+1)
  d <- readSTRef dRef
  return (d,LA.subMatrix (0,0) (d,d) mx)

-- D = P <> A <> Q
-- | Return (d,ker,im) where
-- d: the diagonal entries of Smith normal form
-- ker: basis for the kernel
-- im: basis for the image
kerImOf :: LA.Matrix LA.Z -> (LA.Vector LA.Z,[LA.Vector LA.Z],[LA.Vector LA.Z])
kerImOf mt =
  let (ul,h,ur) = smithNF mt
      (ull,_,ulr) = smithNF ul
      (rk,dmt) = extractMaxNZDiag h
      ker' = ur LA.¿ [rk..(LA.cols ur-1)]
      (_,kert) = hermiteNF (LA.tr' ker')
      im' = (ulr LA.<> ull LA.<> h) LA.¿ [0..rk-1]
      (_,imt) = hermiteNF (LA.tr' im')
  in (LA.takeDiag dmt, LA.toRows kert, LA.toRows imt)

max_stars_repo_name: Junology/linkedgram

max_issues_repo_name: Junology/linkedgram

max_forks_repo_name: Junology/linkedgram


text:
module Main where

import Math.Probably.Sampler
import Math.Probably.RandIO
import Math.Probably.FoldingStats
import qualified Math.Probably.PDF as PDF

import Numeric.LinearAlgebra

import Math.Probably.IterLap

import Graphics.Gnewplot.Exec
import Graphics.Gnewplot.Types
import Graphics.Gnewplot.Style
import Graphics.Gnewplot.Panels
import Graphics.Gnewplot.Instances
import Graphics.Gnewplot.Histogram
import Math.Probably.MCMC

--cov = (2><2) [0.1,0.08,0.08,0.1]

--mn = (2 |> [1,1])

cov3 = (3><3) [0.1,  0.02, 0.08,
               0.02, 0.1,  0.05,
               0.08, 0.05, 0.1]

mn3' = (3 |> [1,1,1])

mn3 = (3 |> [2,2,2])

main = runRIO $ do
  (mn', cov') <- sample $ iterLap [100, 100, 100] (PDF.multiNormal mn3 cov3) (mn3', cov3)
  io $ print mn'
  io $ print cov'

{-  iniampar <- sample $ initialAdaMetWithCov 500 (PDF.multiNormal mn3 cov3) cov' mn3'
  vsamples <- runAdaMetRIO 1500 True iniampar $ (PDF.multiNormal mn3 cov3)
  let (mean,sd) =  (both meanF stdDevF) `runStat` vsamples
  io $ print mean
  io $ print sd -}

--  io $ print ws

max_stars_repo_name: glutamate/probably-base

max_issues_repo_name: glutamate/probably

max_forks_repo_name: glutamate/probably