text:
{-# LANGUAGE BangPatterns #-}
module General.UserNumber
(GeneralSimpleNumber(..)
, ResolveException(..)
, ParseException(..)
, GRatio(..)
, sGR
, reduce
, GeneralRealNumber(..)
, GeneralComplex(..)
, GeneralNumber(..)
, PossNumber(..)
, numShow
, genExp
, realMod
, genRealtoFloat
, gReal
, imagUnit
, getIntegral
, getIntGen
, getFloat
, getRatio
) where

import qualified Data.Complex                     as DC
import           Data.Ratio
import qualified Data.Text.Lazy                   as T
import           Math.NumberTheory.Powers.Squares

data ResolveException = ResolveException T.Text
    deriving Show

data ParseException   = ParseException T.Text Integer
    deriving Show

-- Complex constructor: a:+ b, i.e. 1 + 1j = 1 :+ 1
-- Fraction constructor a % b, i.e. 1 / 2 = 1 % 2 = 2 % 4

--General Number = Either (Either Integer Fractional) (Either Float )

data GeneralSimpleNumber = GInt !Integer | GFlo !Float
  deriving (Show, Eq)

--Show for when printing out. Left as-is for debugging
--instance Show GeneralSimpleNumber where
--  show GInt n = show n
--  show GFlo n = show n

instance Num GeneralSimpleNumber where
  a + b
    | GInt x <- a, GInt y <- b = GInt $ x + y
    | GInt x <- a, GFlo y <- b = GFlo $ fromIntegral x + y
    | GFlo x <- a, GInt y <- b = GFlo $ x + fromIntegral y
    | GFlo x <- a, GFlo y <- b = GFlo $ x + y

  a * b
    | GInt x <- a, GInt y <- b = GInt $ x * y
    | GInt x <- a, GFlo y <- b = GFlo $ fromIntegral x * y
    | GFlo x <- a, GInt y <- b = GFlo $ x * fromIntegral y
    | GFlo x <- a, GFlo y <- b = GFlo $ x * y

  negate (GInt a) = GInt (negate a)
  negate (GFlo a) = GFlo (negate a)

  abs (GInt a) = GInt (abs a)
  abs (GFlo a) = GFlo (abs a)

  signum a
    | GInt x <- a = GInt $ signum x
    | GFlo x <- a, signum x < 0.0 = GInt (-1)
    | GFlo x <- a, signum x > 0.0 = GInt 1
    | GFlo _ <- a, True           = GInt 0

  fromInteger = GInt

instance Ord GeneralSimpleNumber where
  a <= b
    | GInt x <- a, GInt y <- b = x <= y
    | GInt x <- a, GFlo y <- b = fromIntegral x <= y
    | GFlo x <- a, GInt y <- b = x <= fromIntegral y
    | GFlo x <- a, GFlo y <- b = x <= y

data GRatio = GR !Integer !Integer deriving (Show)

instance Eq GRatio where
  GR a b == GR c d = a * d == b * c

instance Num GRatio where
  GR a b + GR c d = ratReduce (a*d+b*c) (b*d)

  GR a b * GR c d = ratReduce (a*b) (c*d)

  negate (GR a b) = GR (negate a) b

  abs (GR a b) = GR (abs a) b

  signum (GR a _) = GR (signum a) 1

  fromInteger n = GR (fromInteger n) 1

instance Fractional GRatio where
  recip (GR a b) = GR b a

  fromRational n = GR (numerator n) (denominator n)

instance Real GRatio where
  toRational (GR a b) = a % b

instance RealFrac GRatio where
  properFraction (GR a b) = (fromInteger $ a `div` b, GR (a `mod` b) b)

instance Ord GRatio where
  GR a b <= GR c d = a*d <= b*c

sGR :: Integer -> Integer -> GRatio
sGR a b = ratReduce (a * signum b) (abs b)

ratReduce :: Integer -> Integer -> GRatio
ratReduce !a !b = let d = gcd a b
               in GR (quot a d) (quot b d)

data GeneralRealNumber = GSimp GeneralSimpleNumber | GRat GRatio deriving (Show, Eq)

reduce :: GRatio -> GeneralRealNumber
reduce r
  | GR a 1 <- r = GSimp (GInt a)
  | otherwise = GRat r

instance Num GeneralRealNumber where
  a + b
    | GSimp x <- a, GSimp y <- b = GSimp $ x + y
    | GSimp (GInt x) <- a, GRat (GR y z) <- b = reduce $ ratReduce (x*z+y) z
    | GSimp (GFlo x) <- a, GRat (GR y z) <- b = GSimp $ GFlo $ x+fromIntegral y/fromIntegral z
    | GRat x <- a, GRat y <- b = reduce $ x + y
    | otherwise = b + a

  a * b
    | GSimp x <- a, GSimp y <- b = GSimp $ x * y
    | GSimp (GInt x) <- a, GRat (GR y z) <- b = reduce $ ratReduce (x*y) z
    | GSimp (GFlo x) <- a, GRat (GR y z) <- b = GSimp $ GFlo $ x*fromIntegral y/fromIntegral z
    | GRat x <- a, GRat y <- b = reduce $ x * y
    | otherwise = b * a

  negate (GSimp a) = GSimp (negate a)
  negate (GRat a)  = GRat (negate a)

  abs (GSimp a) = GSimp (abs a)
  abs (GRat a)  = GRat (abs a)

  signum (GSimp a)       = GSimp (signum a)
  signum (GRat (GR a _)) = GSimp $ GInt $ signum a

  fromInteger n = GSimp (fromInteger n)

instance Fractional GeneralRealNumber where
  recip a
    |GSimp (GInt x) <- a = GRat $ sGR 1 x
    |GSimp (GFlo x) <- a = GSimp $ GFlo $ recip x
    |GRat x <- a = reduce $ recip x

  fromRational n = GRat $ sGR (numerator n) (denominator n)

instance Real GeneralRealNumber where
  toRational (GSimp (GInt x)) = x % 1
  toRational (GSimp (GFlo x)) = toRational x
  toRational (GRat r)         = toRational r

fracHelper :: (RealFrac a, Integral b) => a -> (a -> c) -> (b, c)
fracHelper inner fun = (int, fun frac)
  where
    (int, frac) = properFraction inner

instance RealFrac GeneralRealNumber where
  properFraction (GSimp (GInt x)) = (fromInteger x, GSimp $ GInt 0)
  properFraction (GSimp (GFlo x)) = fracHelper x (GSimp . GFlo)
  properFraction (GRat r)         = fracHelper r GRat

instance Ord GeneralRealNumber where
  GSimp a <= GSimp b = a <= b
  GSimp (GInt a) <= GRat (GR b c)  = a*c <= b
  GSimp (GFlo a) <= GRat (GR b c)  = a * fromInteger c <= fromInteger b
  GRat (GR b c)  <= GSimp (GInt a) = b <= a*c
  GRat (GR b c)  <= GSimp (GFlo a) = fromInteger b <= a * fromInteger c
  GRat a <= GRat b = a <= b

intPow :: GeneralRealNumber -> Integer -> GeneralRealNumber
intPow a n
  | n < 0 = intPow (recip a) (negate n)
  | n == 0 = GSimp $ GInt 1
  | otherwise = a * intPow a (n-1)

expAlgR :: GeneralRealNumber -> Integer -> GeneralRealNumber
expAlgR _ 0 = GSimp $ GInt 1
expAlgR x 1 = x
expAlgR x n
  |n < 0 = expAlgR (recip x) (negate n)
  |even n = expAlgR (x*x) (n `div` 2)
  |otherwise = x * expAlgR (x*x) (n `div` 2)

realExp :: GeneralRealNumber -> GeneralRealNumber -> GeneralRealNumber
realExp _ (GSimp (GInt 0)) = GSimp $ GInt 1
realExp a@(GSimp (GInt _)) (GSimp (GInt y)) = expAlgR a y
realExp (GSimp (GInt x)) (GSimp (GFlo y)) = GSimp $ GFlo $ fromIntegral x ** y
realExp (GSimp (GFlo x)) (GSimp (GInt y)) = GSimp $ GFlo $ x ** fromIntegral y
realExp (GSimp (GFlo x)) (GSimp (GFlo y)) = GSimp $ GFlo $ x**y
realExp r@(GRat (GR _ _)) (GSimp (GInt n)) = expAlgR r n
realExp (GRat (GR a b)) (GSimp c) = realExp (GSimp $ GFlo $ fromInteger a/fromInteger b) (GSimp c)
realExp (GSimp c) (GRat (GR a b)) = realExp (GSimp c) (GSimp $ GFlo $ fromInteger a/fromInteger b)
realExp (GRat (GR a b)) (GRat (GR c d)) = realExp (GSimp $ GFlo $ fromInteger a/fromInteger b) (GSimp $ GFlo $ fromInteger c/fromInteger d)

generalRSqrt :: GeneralRealNumber -> GeneralRealNumber
generalRSqrt (GSimp (GInt a))
  |isSquare a = GSimp $ GInt $ integerSquareRoot a
  |otherwise = GSimp $ GFlo $ sqrt $ fromIntegral a
generalRSqrt (GSimp (GFlo a)) = GSimp $ GFlo $ sqrt a
generalRSqrt (GRat (GR a b))
  |isSquare a && isSquare b = GRat $ GR (integerSquareRoot a) (integerSquareRoot b)
  |otherwise = GSimp $ GFlo $ sqrt $ fromIntegral a/fromIntegral b

genRealtoFloat :: GeneralRealNumber -> Float
genRealtoFloat (GSimp (GInt a)) = fromInteger a
genRealtoFloat (GSimp (GFlo a)) = a
genRealtoFloat (GRat (GR a b))  = fromInteger a/fromInteger b

data GeneralComplex = GC !GeneralRealNumber !GeneralRealNumber
  deriving (Show)

instance Eq GeneralComplex where
  GC a b == GC c d = (a==c) && (b==d)

realPart :: GeneralComplex -> GeneralRealNumber
realPart (GC a _) = a

imagPart :: GeneralComplex -> GeneralRealNumber
imagPart (GC _ b) = b

instance Num GeneralComplex where
  (GC a b) + (GC c d) = GC (a+c) (b+d)

  GC a b * GC c d = GC (a*c-b*d) (a*d+b*c)

  negate (GC a b) = GC (negate a) (negate b)

  abs (GC a b) = GC (generalRSqrt ((a `intPow` 2)+(b `intPow` 2))) (GSimp (GInt 0))

  signum (GC a b) = GC (a / mag) (b / mag)
    where
      mag = realPart $ abs $ GC a b

  fromInteger n = GC (GSimp (GInt n)) 0

instance Fractional GeneralComplex where
  recip (GC a b) = GC (a/((a `intPow` 2)+(b `intPow` 2))) (-b/((a `intPow` 2)+(b `intPow` 2)))

  fromRational n = GC (reduce (ratReduce (numerator n) (denominator n))) (GSimp (GInt 0))

data GeneralNumber = GReal GeneralRealNumber | GComp GeneralComplex
  deriving (Eq)

simpComplex :: GeneralComplex -> GeneralNumber
simpComplex c@(GC a b)
    | b == 0 || b == 0.0 = GReal a
    | otherwise = GComp c

instance Num GeneralNumber where
  a + b
    | GReal x <- a, GReal y <- b = GReal $ x + y
    | GReal x <- a, GComp (GC y z) <- b = simpComplex $ GC (x+y) z
    | GComp x <- a, GComp y <- b = simpComplex $ x + y
    | otherwise = b + a

  a * b
    | GReal x <- a, GReal y <- b = GReal $ x * y
    | GReal x <- a, GComp (GC y z) <- b = simpComplex $ GC (x*y) $ x*z
    | GComp x <- a, GComp y <- b = simpComplex $ x * y
    | otherwise = b * a

  negate (GReal a) = GReal (negate a)
  negate (GComp a) = GComp (negate a)

  abs (GReal a) = GReal (abs a)
  abs (GComp a) = simpComplex (abs a)

  signum (GReal a) = GReal (signum a)
  signum (GComp a) = GComp (signum a)

  fromInteger n = GReal $ GSimp $ GInt n

expAlgC :: GeneralComplex -> Integer -> GeneralNumber
expAlgC _ 0 = GReal $ GSimp $ GInt 1
expAlgC x 1 = GComp x
expAlgC x n
  |n < 0 = expAlgC (recip x) (negate n)
  |even n = expAlgC (x*x) (n `div` 2)
  |otherwise = GComp x * expAlgC (x*x) (n `div` 2)

instance Fractional GeneralNumber where
  recip a
    |GReal x <- a = GReal $ recip x
    |GComp x <- a = GComp $ recip x

  fromRational n = GReal $ fromRational n

genExp :: GeneralNumber -> GeneralNumber -> GeneralNumber
genExp (GReal n) (GReal m) = GReal $ realExp n m
genExp (GComp a) (GReal (GSimp (GInt n))) = expAlgC a n
genExp (GComp a) (GReal (GSimp (GFlo flo))) = GComp (GC (GSimp $ GFlo $ DC.realPart builtInPow) (GSimp $ GFlo $ DC.imagPart builtInPow))
  where
    buildInVer = genRealtoFloat (realPart a) DC.:+ genRealtoFloat (imagPart a)
    builtInPow = buildInVer ** (flo DC.:+ 0)
genExp a@(GComp _) (GReal r@(GRat _)) = genExp a (GReal $ GSimp $ GFlo $ genRealtoFloat r)
genExp (GComp (GC a b)) (GComp (GC c d)) = GComp $ GC (GSimp $ GFlo $ DC.realPart passed) $ GSimp $ GFlo $ DC.imagPart passed
  where
    passed = (genRealtoFloat a DC.:+ genRealtoFloat b) ** (genRealtoFloat c DC.:+ genRealtoFloat d)
genExp (GReal n) j = genExp (GComp (GC n 0)) j

realMod :: GeneralRealNumber -> GeneralRealNumber -> GeneralRealNumber
realMod a (GSimp (GInt 0)) = a
realMod a b = a - (b * f)
  where
    realfloor :: GeneralRealNumber -> Integer
    realfloor (GSimp (GInt x)) = x
    realfloor (GSimp (GFlo x)) = floor x
    realfloor (GRat (GR x y))  = floor $ genRealtoFloat $ GRat $ GR x y
    f = GSimp $ GInt $ realfloor $ a/b

realShow :: GeneralRealNumber -> String
realShow num
  |GRat (GR x y) <- num = show x ++ "/" ++ show y
  |GSimp (GInt x) <- num = show x
  |GSimp (GFlo x) <- num = show x

numShow :: GeneralNumber -> String
numShow num
  |GComp (GC (GSimp (GInt 0)) y) <- num = realShow y ++ "j"
  |GComp (GC x y) <- num = realShow x ++ "+" ++ realShow y ++ "j"
  |GReal x <- num = realShow x

instance Show GeneralNumber where
  show = numShow

data PossNumber = Kept GeneralNumber
                | Dropped GeneralNumber
  deriving (Eq)

instance Show PossNumber where
  show (Kept n)    = show n
  show (Dropped n) = "//" ++ show n ++ "//"

gReal :: (Integral a) => a -> GeneralNumber
gReal n = GReal $ GSimp $ GInt $ fromIntegral n

imagUnit = GComp $ GC (GSimp $ GInt 0) (GSimp $ GInt 1)

getIntegral :: (Integral a) => GeneralNumber -> Maybe a
getIntegral (GReal (GSimp (GInt n))) = Just $ fromInteger n
getIntegral _                        = Nothing

getIntGen :: (Integer -> a) -> GeneralNumber -> Maybe a
getIntGen p (GReal (GSimp (GInt n))) = Just $ p n
getIntGen _ _                        = Nothing

getFloat :: GeneralNumber -> Maybe Float
getFloat (GReal (GSimp (GFlo f))) = Just f
getFloat _                        = Nothing

getRatio :: GeneralNumber -> Maybe Rational
getRatio (GReal (GRat (GR p q))) = Just $ p % q
getRatio _                       = Nothing

max_stars_repo_name: cwstra/rpdiscordrewrite

max_issues_repo_name: cwstra/rpdiscordrewrite

max_forks_repo_name: cwstra/rpdiscordrewrite


text:
import Data.Map.Strict ((!))
import Numeric.LinearAlgebra hiding ((!))
import LeastSquares

leastSquaresEstimate :: Vector R
leastSquaresEstimate = result ! "x"
  where result = minimize (SumSquares x) [c * x :==: d]
        c = Mat $ (2><3) [3.0, 7.0, -1.5, 2.8, 10.0, -5.3]
        d = Vec $ vector [-1.0, 3.5]
        x = Var "x" 3

ext: hs

max_stars_repo_name: tridao/leastsquares

max_issues_repo_name: tridao/leastsquares

max_forks_repo_name: tridao/leastsquares


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE TemplateHaskell   #-}
{-# LANGUAGE UnicodeSyntax     #-}
{-# LANGUAGE BangPatterns      #-}
module NLP.Albemarle.LSA
    (-- * Creating a model
    lsa
    -- * Managing terms
    , rebase
    -- * Managing topics
    , pad
    , trim
    -- * Internal methods (may be subject to change)
    , SVD.sparsify
    , SVD.sparseSvd
    , batchLSA
    -- * Types (may be subject to change)
    , LSAModel(..)
    , termvectors
    , topicweights
    ) where
import NLP.Albemarle
import NLP.Albemarle.Dict (Dict)
import qualified NLP.Albemarle.Dict as Dict
import ClassyPrelude hiding (Vector)
import qualified Data.Vector.Generic as Vec
import qualified Data.Vector.Unboxed as UVec
import qualified Data.Vector.Storable as SVec
import qualified Data.HashMap.Strict as HashMap
import Numeric.LinearAlgebra (Vector, Matrix, tr, diag, (|||), (===))
import qualified Numeric.LinearAlgebra as HMatrix
import qualified Numeric.LinearAlgebra.Devel as HMatrix
import qualified Numeric.LinearAlgebra.SVD.SVDLIBC as SVD
import qualified System.IO.Streams as Streams
import System.IO.Streams (Generator, InputStream, OutputStream)
import Lens.Micro.TH
import Lens.Micro
import Data.Tuple

-- | An LSA model (The singular values and right singular vectors of truncated
--   SVD on a term-document matrix, where documents are rows and terms columns)
data LSAModel = LSAModel {
  _topicweights :: !(HMatrix.Vector Double), -- ^ Topic weights
  _termvectors :: !(Matrix Double) -- ^ Rows are topics, columns are terms
} deriving (Show, Eq)
makeLenses ''LSAModel

-- This instance may be a fib. I'm not sure this is actually associative.
-- Even if it is, I imagine there are practical numerical stability concerns.
-- TODO: In particular, I think we need to weight the left and weight according
-- to how many documents they represent.
instance Monoid LSAModel where
  mempty = LSAModel mempty mempty
  mappend left right
    | left == mempty = right
    | right == mempty = left
    | otherwise = trim target_len $ pad target_len $ LSAModel s v
    where
      target_len = max (topicCount left) (topicCount right)
      (s, v) = HMatrix.rightSV $ combine left === combine right
instance Semigroup LSAModel where
  (<>) = mappend

-- | Truncate an LSAModel to a specific number of topics, if there are too many
trim :: Int -- ^ The maximum number of topics the model should have (inclusive)
     -> LSAModel -> LSAModel
trim count model
  | topicCount model <= count = model
  | otherwise = LSAModel
    (Vec.take count $ model^.topicweights)
    (HMatrix.takeRows count $ model^.termvectors)

-- | Pad an LSAModel with extra empty topics of there are too few
pad :: Int -- ^ The minimum number of topics the model should have (inclusive)
    -> LSAModel -> LSAModel
pad count model
  | topicCount model >= count = model
  | otherwise = let
    (height, width) = HMatrix.size $ model^.termvectors
    missing = count - height
    in LSAModel
      (model^.topicweights <> Vec.replicate missing 0)
      (model^.termvectors === HMatrix.konst 0 (missing, width))

-- | Multiply the vectors and values of an SVD to make one matrix
combine :: LSAModel -> Matrix Double
combine model = diag (model^.topicweights) <> (model^.termvectors)

sparseToCSR :: SparseMatrix -> HMatrix.CSR
sparseToCSR (SparseMatrix width vecs) =
  let
    docs = (Vec.fromList . (\ (SparseVector len wds) -> wds) <$> vecs)
    convI = Vec.map fromIntegral . Vec.convert
    counts = Vec.map snd <$> docs         :: [UVec.Vector Double]
    concatcounts = Vec.concat counts      ::  UVec.Vector Double
    storcounts = Vec.convert concatcounts ::  SVec.Vector Double
  in HMatrix.CSR {
    HMatrix.csrVals = storcounts,
    HMatrix.csrCols = convI
      $ Vec.concat
      $ Vec.map ((+1).fst)
      <$> docs,
    HMatrix.csrRows = convI
      $ Vec.scanl (+) 1
      $ Vec.fromList
      $ Vec.length
      <$> docs,
    HMatrix.csrNCols = width,
    HMatrix.csrNRows = length docs
  }

-- | Generate an LSA Model with N topics, based on a sparse matrix.
--
--   You can make the sparse matrix with a Dict and some documents.
--   See Dict.asSparseMatrix.
lsa :: Int -- ^ Number of vectors/dimensions in the new space
    -> SparseMatrix -- ^ Sparse representation of the term-document matrix
    -> LSAModel -- ^ A dense representation of the term-topic matrix
lsa top_vectors termdoc =
  -- This astounds me, but the number of singular values may not match the
  -- number of singular vectors. We have to balance before padding.
  pad top_vectors $ case compare missing_values 0 of
    LT -> LSAModel (Vec.take top_vectors s) vt
    EQ -> LSAModel s vt
    GT -> LSAModel (s <> Vec.replicate missing_values 0) vt
  where
    (u, s, vt) = SVD.sparseSvd top_vectors $ sparseToCSR termdoc
    missing_values = fst (HMatrix.size vt) - Vec.length s

-- | Rebase an LSA model made with one Dict to one made with another Dict.
--
--   By using this, you can split a corpus into batches, make Dicts and topic
--   models, then merge the Dicts, rebase the Models to the merged Dict, and
--   then merge the models. It allows working in parallel, and limits memory
--   usage.
rebase :: Dict -> Dict -> LSAModel -> LSAModel
rebase d1 d2 model
  | model == mempty = mempty
  | otherwise = LSAModel {
    _termvectors = (HMatrix.¿) (model^.termvectors) $ Dict.select d1 d2,
    _topicweights = model^.topicweights
  }

-- | Get the number of topics in a model
topicCount :: LSAModel -> Int
topicCount model = Vec.length $ model^.topicweights

-- | SVD with some transposes, for convenience and speed
-- Normally (U, Sigma, V^T) = svd A, but then the rows are topics and the cols
-- are the documents/words (in U and V). But we use C-style (row major) matrices
-- which means for large matrices, getting the vector of one word or one doc
-- will require reading the whole model from memory. (A terrible waste of cache)
-- So instead we use (U^T, sigma, V) where the rows are the vector embeddings
-- of documents and words instead.
--
-- Prefer lsa (rather than batchLSA) - batchLSA is experimental and may be
-- removed
batchLSA :: Int -> HMatrix.CSR -> (Matrix Double, HMatrix.Vector Double, Matrix Double)
batchLSA top_vectors csr = (tr u, s, tr vt)
  where (u, s, vt) = SVD.sparseSvd top_vectors csr

max_stars_repo_name: SeanTater/albemarle

max_issues_repo_name: SeanTater/albemarle

max_forks_repo_name: SeanTater/albemarle


text:
{-|
Module      : Mandelbrot
Description : Algorithm of Mandelbrot functions

Module which get resolution of image and points and calculate
correct color to point based on Mandelbrot algorithm
-}
module Mandelbrot (drawMandelbrot, runMandelbrotTests, runQuickCheckMandelbrot) where
import Data.Complex
import Graphics.GD
import Plot
import Test.HUnit
import Test.QuickCheck

-- |Max number of iteration to approximate set
maxIteration = 150

-- |Main function, invoke draw Mandelbrot set
drawMandelbrot :: Size -- ^ Resolution of generating image
                -> Bound -- ^ Four coordinates of funtion bound
                -> String -- ^ Name of generating file
                -> IO() -- ^ Result: Save file
drawMandelbrot size bound fileName = drawPlot drawMandelbrot_ size bound fileName

-- |Generating colors for points
drawMandelbrot_ :: Coord -- ^ Point to draw
                -> Color -- ^ Color of the point calculated from mandelbrot formula
drawMandelbrot_ (x,y) = rgb (floor $ f * 255) (floor $ f * 120) 0
                        where i = mandelbrot (x :+ y)
                              f = fromIntegral i / fromIntegral maxIteration

-- |Recursive formula for points
mandelbrot :: Complex Double -- ^ Checking point
              -> Int -- ^ Number of iteration
mandelbrot z = length $ takeWhile (\c -> magnitude c <= 2) $ take maxIteration $ iterate (\w -> w^2 + z) 0

test1 = TestCase (assertEqual "color1" 16742400 (drawMandelbrot_ (0,0)))
test2 = TestCase (assertEqual "color2" 196864 (drawMandelbrot_ (1,1)))
test3 = TestCase (assertEqual "color3" 65536 (drawMandelbrot_ (10,10)))
test4 = TestCase (assertEqual "mandelbrot1" 2 (mandelbrot (1 :+1)))
test5 = TestCase (assertEqual "mandelbrot2" 1 (mandelbrot (0 :+ 20)))
test6 = TestCase (assertEqual "mandelbrot3" 1 (mandelbrot (123 :+ 0)))
test7 = TestCase (assertEqual "mandelbrot3" 150 (mandelbrot (0 :+ 0)))

tests = TestList [TestLabel "test1" test1, TestLabel "test2" test2, TestLabel "test3" test3, TestLabel "test4" test4, TestLabel "test5" test5, TestLabel "test6" test6, TestLabel "test7" test7]

-- |Run tests
runMandelbrotTests :: IO Counts
runMandelbrotTests = do runTestTT tests

max_Color :: Coord -> Bool
max_Color (x,y) =   drawMandelbrot_ (x,y)  <= 16777216

min_Color :: Coord -> Bool
min_Color (x,y) =   drawMandelbrot_ (x,y)  >= 0

-- |Run quick tests
runQuickCheckMandelbrot = do quickCheck max_Color >> quickCheck min_Color

max_stars_repo_name: jsamol/fractal-painter

max_issues_repo_name: jsamol/fractal-painter

max_forks_repo_name: jsamol/fractal-painter


text:
{-- --:= Mark: find better way to construct complex number
type ImaginaryUnit = Integer
type RealNumber = Integer
data ComplexNumber = ComplexN {real :: RealNumber,
                               imag :: ImaginaryUnit}
--}

import           Data.Complex

readComplex :: String -> Complex Double
readComplex = read

result :: (String, String) -> Complex Double
result (x,y) = (readComplex x) * (readComplex y)

main = do
  let exm1 = ("1 :+ 1","1 :+ 1")
  let exm2 = ("1 :+ -1","1 :+ -1")
  print $ result exm1
  print $ result exm2

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
{-# LANGUAGE OverlappingInstances #-}
{-# LANGUAGE UnicodeSyntax #-}

-- thanks to http://www.muitovar.com/gtk2hs/app1.html

--module Test where

import Control.Concurrent
import Control.Concurrent.MVar

import Control.Monad.Trans

import Graphics.UI.Gtk hiding(Circle,Cross)
import qualified Graphics.Rendering.Cairo as C
import qualified Graphics.Rendering.Pango as P

import Data.Colour.Names

import Data.Packed.Vector
--import Data.Packed.Random
import Data.Packed()

--import Prelude.Unicode

import qualified Data.Array.IArray as A

import Numeric.LinearAlgebra

import Numeric.GSL.Statistics

import Graphics.Rendering.Plot

import Debug.Trace

ln = 25
ts = linspace ln (0,1)
rs = ln |> take ln [0.306399512330476,-0.4243863460546792,-0.20454667402138094,-0.42873761654774106,1.3054721019673694,0.6474765138733175,1.1942346875362946,-1.7404737823144103,0.2607101951530985,-0.26782584645524893,-0.31403631431884504,3.365508546473985e-2,0.6147856889630383,-1.191723225061435,-1.9933460981205509,0.6015225906539229,0.6394073044477114,-0.6030919788928317,0.1832742199706381,0.35532918011648473,0.1982646055874545,1.7928383756822786,-9.992760294442601e-2,-1.401166614128362,-1.1088031929569364,-0.827319908453775,1.0406363628775428,-0.3070345979284644,0.6781735212645198,-0.8431706723519456,-0.4245730055085966,-0.6503687925251668,-1.4775567962221399,0.5587634921497298,-0.6481020127107823,7.313441602898768e-2,0.573580543636529,-0.9036472376122673,2.650805059813826,9.329324044673039e-2,1.9133487025468563,-1.5366337588254542,-1.0159359710920388,7.95982933517428e-2,0.5813673663649735,-6.93329631989878e-2,1.1024137719307867,-0.6046286796589855,-0.8812842030098401,1.4612246471009083,0.9584060744500491,9.210899579679932e-2,-0.15850413664405813,-0.4754694827227343,0.8669922262489788,0.4593351854708853,-0.2015350278936992,0.8829710664887649,0.7195048491420026]

ss = sin (15*2*pi*ts)
ds = 0.25*rs + ss
es = constant (0.25*(stddev rs)) ln

fs :: Double -> Double
fs = sin . (15*2*pi*)

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

pts = linspace 1000 (0 :: Double,10*pi)
fx = (\t -> t * sin t) pts
fy = (\t -> t * cos t) pts

hx = fromList [1,3,5,8,11,20,22,26] :: Vector Double
hy = fromList [10,11,15,17,14,12,9] :: Vector Double

lx = fromList [1,2,3,4,5,6,7,8,9,10] ∷ Vector Double
ly = fromList [50000,10000,5000,1000,500,100,50,10,1] ∷ Vector Double

mx = linspace 100 (1,10) ∷ Vector Double
my = linspace 100 (1,10000) ∷ Vector Double

cx = fromList [1,2,3,4,5] ∷ Vector Double
cyl = fromList [8,10,12,13,8] ∷ Vector Double
cyu = fromList [10,12,16,5,10] ∷ Vector Double
cel = cyl - 1
ceu = cyu + 1

at = linspace 1000 (0,2*pi) ∷ Vector Double
ax = sin at

figure = do
--         setPlots 1 1
{-
         withPlot (1,1) $ do
                          setDataset [(Hist,hx,hy)]
                          addAxis XAxis (Side Lower) $ return ()
                          addAxis YAxis (Side Lower) $ return ()
-}{-                          setRange XAxis Lower (-4*pi) (1*pi)
                          setRange YAxis Lower (-4*pi) (1*pi) -}
{-                          setRange XAxis Lower 0 32
                          setRange YAxis Lower 0 20
-}
        withLineDefaults $ setLineWidth 2
        withTextDefaults $ setFontFamily "OpenSymbol"
        withTitle $ setText "Testing plot package:"
        withSubTitle $ do
                       setText "with 1 second of a 15Hz sine wave"
                       setFontSize 10
        setPlots 1 1

        withPlot (1,1) $ do

--                         setDataset (ts,[area ds blue])
--                         setDataset (ts,[impulse fs blue])
--                         setDataset (ts,[point (ds,es,"Sampled data") (Bullet,green)
--                                        ,line (fs,"15 Hz sinusoid") blue])
--                         setDataset [(Line,fx,fy)]
--                         setDataset (ts,[bar (ds,"Sampled data") (10 :: Double,green,3:: Double,blue)
--                                        ,line (fs,"15 Hz sinusoid") blue])
--                         setDataset [(Line,mx,my)]
--                         setDataset (Whisker,cx,[((cyl,cyu),(cel,ceu))])
                         setDataset (Line,at,[ax])
                         addAxis XAxis (Side Lower) $ do
                                                      setGridlines Major True
                                                      withAxisLabel $ setText "time (s)"
                         addAxis YAxis (Side Lower) $ do
                                                      setGridlines Major True
                                                      withAxisLabel $ setText "amplitude (α)"
--                         addAxis XAxis (Value 0) $ return ()
                         setRangeFromData XAxis Lower Linear
                         setRangeFromData YAxis Lower Linear
                         withAnnotations $ do
                           arrow True (pi/2,0.5) (0,0) (return ())
                           oval True (1.5,0) (pi,0.5) $ setBarColour blue
                           rect True (0.5,0.5) (2,0.6) $ (return ())
                           glyph (4,0.2) (return ())
                           text (3,0.0) (setText "from the α to the Ω")
                           cairo (\_ _ _ _ -> do
                                    C.newPath
                                    C.moveTo 3 0.75
                                    C.lineTo 4 (-0.5)
                                    C.stroke
                                    C.rectangle (pi/2) (-1) (pi/2) 1
                                    C.fill
                                 )
--                         setRange YAxis Lower Log (-1.25) 1.25
--                         setLegend True NorthEast Inside
--                         withLegendFormat $ setFontSize 6
{-
         withPlot (1,1) $ do
                          setDataset (ident 300 :: Matrix Double) --ms
                          addAxis XAxis (Side Lower) $ setTickLabelFormat "%.0f"
                          addAxis YAxis (Side Lower) $ setTickLabelFormat "%.0f"
                          setRangeFromData XAxis Lower
                          setRangeFromData YAxis Lower
-}

display :: ((Int,Int) -> C.Render ()) -> IO ()
display r = do
   initGUI       -- is start

   window <- windowNew
   set window [ windowTitle := "Cairo test window"
              , windowDefaultWidth := 600
              , windowDefaultHeight := 400
              , containerBorderWidth := 1
              ]

--   canvas <- pixbufNew ColorspaceRgb True 8 300 200
--   containerAdd window canvas
   frame <- frameNew
   containerAdd window frame
   canvas <- drawingAreaNew
   containerAdd frame canvas
   widgetModifyBg canvas StateNormal (Color 65535 65535 65535)

   widgetShowAll window

   on canvas exposeEvent $ tryEvent $ do s <- liftIO $ widgetGetSize canvas
                                         drw <- liftIO $ widgetGetDrawWindow canvas
                                         --dat <- liftIO $ takeMVar d
                                         --liftIO $ renderWithDrawable drw (circle 50 10)
                                         liftIO $ renderWithDrawable drw (r s)

   onDestroy window mainQuit
   mainGUI

main = display $ render figure

test = writeFigure PDF "test.pdf" (400,400) figure

max_stars_repo_name: tsurucapital/plot

max_issues_repo_name: tsurucapital/plot

max_forks_repo_name: tsurucapital/plot


text:
module Main where
import Plot
import Graphics.GD
import Data.Complex

imageDim :: Size
imageDim = (600,600)

-- How many iterations before  we declare the point doesn't diverge
maxIter = 127

-- Coordinates in the complex plane to graph
-- First coordinate must be top left.
mbWindow :: Window
mbWindow = ( (-1.9, 1.3), (0.7, -1.3) ) -- Full view

main :: IO ()
main = drawPlot drawMandelbrot imageDim mbWindow "mandelbrot.png"

drawMandelbrot :: Coordinate -> Color
drawMandelbrot (x,y) = colorIterations $ mandelbrot (x :+ y) (0 :+ 0) 0

-- Count number of iterations for the given point to diverge
-- Non-diverging points return 0
mandelbrot :: Complex Double -- Coordinate to test
           -> Complex Double -- Iterating Z value
           -> Int -- Current iteration
           -> Int -- Iterations before diverging
mandelbrot c z iter
    | iter > maxIter = 0
    | otherwise = let z' = z^2 + c in
                  if magnitude z' > 2
                  then iter
                  else mandelbrot c z' (iter+1)

colorIterations :: Int -> Color
colorIterations x
    | x > maxIter = rgb 255 255 255
    | otherwise = let c = x + 30
                  in rgb c x (x^2)

max_stars_repo_name: typeclasses/fractal-sets

max_issues_repo_name: typeclasses/fractal-sets

max_forks_repo_name: typeclasses/fractal-sets


text:
module EKF
  ( KF
  , kf_x
  , kf_p
  , kf_t

  , makeFilter
  , predict
  , update
  ) where

import Prelude hiding ((<>))

import App (Sensor (..))

import Numeric.LinearAlgebra

data KF = KF
  { kf_x :: Vector Double
  , kf_p :: Matrix Double
  , kf_t :: Word
  } deriving Show

makeFilter (Laser px py t) =
  KF { kf_x = 4 |> [px,
                    py,
                    0,
                    0]
     , kf_p = (4><4) [1, 0,    0,    0,
                      0, 1,    0,    0,
                      0, 0, 1000,    0,
                      0, 0,    0, 1000]
     , kf_t = t }

makeFilter (Radar rho phi rho' t) =
  KF { kf_x = 4 |> [rho  * cos phi,
                    rho  * sin phi,
                    rho' * cos phi,
                    rho' * sin phi]
     , kf_p =  (4><4) [1, 0,    0,    0,
                       0, 1,    0,    0,
                       0, 0, 1000,    0,
                       0, 0,    0, 1000]
     , kf_t = t }

predict :: Word -> KF -> KF
predict t kf =
  KF { kf_x = x
     , kf_p = p
     , kf_t = t }
  where
    -- Time inside the filter.
    filterTime =
      kf_t kf

    -- Delta-time. Essentially a timestep.
    dt =
      (fromIntegral (t - filterTime)) / 1000000.0
    dt2 = dt ** 2
    dt3 = dt ** 3
    dt4 = dt ** 4

    -- Acceleration noise tuning parameters
    ax = 3
    ay = 3

    noiseAx = ax ** 2
    noiseAy = ay ** 2

    -- State transition matrix
    f = (4><4) [1, 0, dt,  0,
                0, 1,  0, dt,
                0, 0,  1,  0,
                0, 0,  0,  1]

    -- Noise covariance matrix
    q = (4><4) [-- Row 1
                noiseAx * 0.25 * dt4, 0, noiseAx * 0.5 * dt3, 0,
                -- Row 2
                0, noiseAy * 0.25 * dt4, 0, noiseAy * 0.5 * dt3,
                --  Row 3
                noiseAx *0.5 * dt3, 0, noiseAx * dt2, 0,
                -- Row 4
                0, noiseAy * 0.5 * dt3, 0, noiseAy * dt2]

    -- Noise vector
    v  = 4 |> [ax * 0.5 * dt2,
               ay * 0.5 * dt2,
               ax * dt,
               ay * dt]

    -- Estimated step based on constant velocity (CV) motion model
    x = f #> kf_x kf + v
    p = f <> kf_p kf <> (tr f) + q

update (Laser px py t) kf =
  KF { kf_x = x'
     , kf_p = p'
     , kf_t = t }
  where
    kf' = predict t kf

    -- Measurment transition matricies
    h = (2><4) [1, 0, 0, 0,
                0, 1, 0, 0]
    r = (2><2) [0.0225, 0.0000,
                0.0000, 0.0225]
    y = 2 |> [px, py] - (h #> kf_x kf')

    -- Kalman Filter Equations
    i  = ident 4
    ht = tr h
    s  = h <> kf_p kf' <> ht + r
    si = inv s
    k  = kf_p kf' <> ht <> si

    -- Estimation based on kalman filter gain
    x' = kf_x kf' + (k #> y)
    p' = (i - k <> h) <> kf_p kf'

update (Radar rho phi rho' t) kf =
  KF { kf_x = x'
     , kf_p = p'
     , kf_t = t }
  where
    kf' = predict t kf

    -- Calculate the Jacobian matrix values
    px = kf_x kf' ! 0
    py = kf_x kf' ! 1
    vx = kf_x kf' ! 2
    vy = kf_x kf' ! 3

    c1 = px ** 2 + py ** 2
    c2 = sqrt c1
    c3 = c1 * c2

    zRho  = c2;
    zPhi  = (atan2 py px)
    zRho' = (px*vx + py*vy) / zRho

    -- Measurment transition matricies
    h = (3><4) [-- Row 1
                px / c2, py / c2, 0, 0,
                -- Row 2
                -py / c1, px / c1, 0, 0,
                -- Row 3
                py*(vx * py - vy * px)/c3,
                px*(vy * px - vx * py)/c3,
                px/c2,
                py/c2]
    r = (3><3) [0.09, 0.0000, 0.00,
                0.00, 0.0009, 0.00,
                0.00, 0.0000, 0.09]
    z = 3 |> [rho, phi, rho'] - 3 |> [zRho, zPhi, zRho']
    y = 3 |> [z ! 0, atan2 (sin (z ! 1)) (cos (z ! 1)), z ! 2]

    -- Kalman filter equations
    i  = ident 4
    ht = tr h
    s  = h <> kf_p kf' <> ht + r
    si = inv s
    k  = kf_p kf' <> ht <> si

    -- Estimation based on kalman filter gain
    x' = kf_x kf' + (k #> y)
    p' = (i - k <> h) <> kf_p kf'


text:
-- ch1 ex. 1.3

import qualified Data.Vector as V

import Statistics.Sample.Histogram (histogram)
import Statistics.Sample (mean, varianceUnbiased, stdDev, range)

-- uranium 258 in soil (picocuries per gram)
samples :: V.Vector Double
samples = V.fromList
  [ 0.74, 6.47, 1.90, 2.69, 0.75
  , 0.32, 9.99, 1.77, 2.41, 1.96
  , 1.66, 0.70, 2.42, 0.54, 3.36
  , 3.59, 0.37, 1.09, 8.32, 4.06
  , 4.55, 0.76, 2.03, 5.70,12.48
  ]

hist :: (V.Vector Double, V.Vector Int)
hist = histogram 8 samples

main :: IO ()
main = do
  let (lowerBs, hs) = hist
  putStrLn $ "The lower bounds for each bin: " ++ show lowerBs
  putStrLn $ "The histogram is: " ++ show hs

  let av = mean samples
      sd = stdDev samples
  putStrLn $ "The average value is: " ++ show av ++ "and the standard deviation is: " ++ show sd

  putStrLn $ "The range/4 is: " ++ show (range samples / 4)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
module Parser where
import Text.ParserCombinators.Parsec
import System.Environment
import Control.Monad
import Utilities.Number
import Numeric
import Data.Ratio
import Data.Complex
import Data.Array
import Utilities.Types

symbol :: Parser Char
symbol = oneOf "!$%&|*+-/:<=>?@^_~"

spaces1 :: Parser ()
spaces1 = skipMany1 space

escapedChars :: Parser Char
escapedChars = do
    char '\\' -- a backslash
    x <- oneOf "\\\"nrt" -- return either backslash or doublequote
    return $ case x of
        '\\' -> x
        '"' -> x
        'n' -> '\n'
        'r' -> '\r'
        't' -> '\t'

parseCharacter :: Parser LispVal
parseCharacter = do
    try $ string "#\\"
    value <- try (string "newline" <|> string "space")
        <|> do { x <- anyChar; notFollowedBy alphaNum; return [x]}
    return $ Character $ case value of
        "space" -> ' '
        "newline" -> '\n'
        _ -> head value

parseBool :: Parser LispVal
parseBool = do
    char '#'
    (char 't' >> return (Bool True)) <|> (char 'f' >> return (Bool False))

parseString :: Parser LispVal
parseString = do
    char '"'
    x <- many $ escapedChars <|> noneOf "\"\\"
    char '"'
    return $ String x

parseAtom :: Parser LispVal
parseAtom = do
    first <- letter <|> symbol
    rest <- many (letter <|> digit <|> symbol)
    let atom = first:rest
    return $ case atom of
        "#t" -> Bool True
        "#f" -> Bool False
        v -> Atom v

parseDecimal1 :: Parser LispVal
parseDecimal1 = Number . read <$> many1 digit

parseDecimal2 :: Parser LispVal
parseDecimal2 = do
    try $ string "#d"
    x <- many1 digit
    (return . Number . read) x

parseHex :: Parser LispVal
parseHex = do
    try $ string "#x"
    x <- many1 hexDigit
    return $ Number (hex2dig x)

parseOct :: Parser LispVal
parseOct = do
    try $ string "#o"
    x <- many1 octDigit
    return $ Number (oct2dig x)

parseBin :: Parser LispVal
parseBin = do
    try $ string "#o"
    x <- many1 (oneOf "10")
    return $ Number (bin2dig x)

parseNumber :: Parser LispVal
parseNumber = parseDecimal1 <|> parseDecimal2 <|> parseHex <|> parseOct <|> parseBin

parseFloat :: Parser LispVal
parseFloat = do
    x <- many1 digit
    char '.'
    y <- many1 digit
    return $ Float (fst . head $ readFloat (x ++ "." ++ y))

parseRatio :: Parser LispVal
parseRatio = do
    x <- many1 digit
    char '/'
    y <- many1 digit
    return $ Ratio (read x % read y)

parseComplex :: Parser LispVal
parseComplex = do
    x <- try parseFloat <|> parseDecimal1 <|> parseDecimal2
    char '+'
    y <- try parseFloat <|> parseDecimal1 <|> parseDecimal2
    char 'i'
    return $ Complex (toDouble x :+ toDouble y)

-- parseList :: Parser LispVal
-- parseList = liftM List $ sepBy parseExpr spaces1

-- parseDottedList :: Parser LispVal
-- parseDottedList = do
--     head <- endBy parseExpr spaces1
--     tail <- char '.' >> spaces1 >> parseExpr
--     return $ DottedList head tail

parseList :: Parser LispVal
parseList = do
    char '(' >> spaces
    head <- parseExpr `sepEndBy` spaces1
    do
        char '.' >> spaces1
        tail <- parseExpr
        spaces >> char ')'
        return $ DottedList head tail
        <|> (spaces >> char ')' >> (return $ List head))

parseQuoted :: Parser LispVal
parseQuoted = do
    char '\''
    x <- parseExpr
    return $ List [Atom "quote", x]

parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
    char '`'
    x <- parseExpr
    return $ List [Atom "quasiquote", x]

parseUnQuote :: Parser LispVal
parseUnQuote = do
    char ','
    x <- parseExpr
    return $ List [Atom "unquote", x]

parseVector :: Parser LispVal
parseVector = do
    arrayValues <- sepBy parseExpr spaces1
    return $ Vector (listArray (0, length arrayValues - 1) arrayValues)

parseExpr :: Parser LispVal
parseExpr = parseAtom
    <|> parseString
    <|> try parseComplex
    <|> try parseFloat
    <|> try parseRatio
    <|> try parseNumber
    <|> parseQuoted
    <|> try (do
        string "#("
        x <- parseVector
        char ')'
        return x)
    <|> parseList
    <|> try parseBool
    <|> try parseCharacter
    <|> parseQuasiQuoted
    <|> parseUnQuote

max_stars_repo_name: jordanosborn/Haskell-Lisp

max_issues_repo_name: jordanosborn/Haskell-Lisp

max_forks_repo_name: jordanosborn/Haskell-Lisp


text:
-- |
-- Module      : Application
-- Description : Utilties for creating GUI applications
-- Copyright   : (c) Jonatan H Sundqvist, year
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created date year

-- TODO | - Refactor (move types to separate module, add lenses)
--        - Input queries, composite events (InputState type?)
--        - Redesign API (less fragile and add-hoc) (divide into app types, eg. full screen canvas, etc.)
--        - Optional event wrappers
--        - Printf-style event coupling (with a Listener class (?))

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC directives
--------------------------------------------------------------------------------------------------------------------------------------------
{-# LANGUAGE RankNTypes #-}

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module Southpaw.Interactive.Application where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Data.Complex --
import Data.Functor --
import Data.IORef   --
import Data.Maybe   --

import Control.Lens hiding (set)
import Control.Monad       (liftM, when, void, forM) --
import Control.Applicative (pure, (<*>))

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

import Southpaw.Interactive.Types
import Southpaw.Interactive.Lenses

--------------------------------------------------------------------------------------------------------------------------------------------
-- Data
--------------------------------------------------------------------------------------------------------------------------------------------

-- |
-- TODO: Rename (eg. no-listeners), create a separate default value (?)
nolisteners :: EventMap self s
nolisteners = EventMap Nothing Nothing Nothing Nothing Nothing Nothing Nothing Nothing

-- |
-- TODO: Rename (eg. no-listeners), create a separate default value (?)
-- defaultListeners :: EventMap s
-- defaultListeners = nolisteners { }

--------------------------------------------------------------------------------------------------------------------------------------------
-- Functions
--------------------------------------------------------------------------------------------------------------------------------------------

-- Master listeners ------------------------------------------------------------------------------------------------------------------------

-- TODO: Hoogle section

-- |
masterdraw :: (s -> Cairo.Render ()) -> IORef s -> Cairo.Render ()
masterdraw draw stateref = do
	state <- Cairo.liftIO $ readIORef stateref
	draw state

-- General utilities -----------------------------------------------------------------------------------------------------------------------

-- |
-- TODO: Make polymorphic
widgetSize :: WidgetClass self => self -> IO (Complex Int)
widgetSize widget = (:+) <$> widgetGetAllocatedWidth widget <*> widgetGetAllocatedHeight widget

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

-- |
-- TODO: Rename (duh)
-- TODO: Refactor
-- TODO: Factor out default event handlers
-- TODO: Use state value from the App (?)
wholecanvas :: WidgetSettings Window state -> WidgetSettings DrawingArea state -> state -> IO (App FullCanvasLayout state)
wholecanvas windowSettings canvasSettings appstate = do
	initGUI

	stateref <- newIORef appstate

	window'  <- makeWidget windowNew stateref $ windowSettings
	frame'   <- frameNew
	canvas'  <- makeWidget drawingAreaNew stateref $ canvasSettings

	containerAdd frame' canvas'
	set window' [ containerChild := frame' ]

	widgetShowAll window'
	return $ App { _window=window', _layout=FullCanvasLayout { _canvas=canvas' }, _appsize=_settingsize windowSettings, _state=stateref }

-- |
-- TODO: Use readstate for everything, or remove it (?)
-- TODO: Find out how to deal with 'onanimate'
-- TODO: Event wrappers (eg. draw)
attachListeners :: WidgetClass self => self -> IORef state -> EventMap self state -> IO ()
attachListeners self stateref eventmap = do
	simpleBind buttonPressEvent   onmousedown
	simpleBind buttonReleaseEvent onmouseup
	simpleBind motionNotifyEvent  onmousemotion
	simpleBind keyPressEvent      onkeydown
	simpleBind keyReleaseEvent    onkeyup
	simpleBind deleteEvent        ondelete
	-- simpleBind draw           (\emap -> maybe Nothing (void . Just . masterdraw dr) (ondraw emap))
	maybe pass (\ondr       -> void $ (self `on` draw $ masterdraw ondr stateref)) (ondraw eventmap)
	maybe pass (\(fps, ona) -> void $ timeoutAdd (ona self stateref) fps)          (onanimate eventmap)

	pass
	-- perhapsBind self key onanimate
	-- onresize, -- configureEvent
	-- ondrawIO
	-- perhapsBind self  onquit
	-- forM [(self `on` motionNotifyEvent,  onmousemotion),
	--       (self `on` buttonPressEvent,   onmousedown),
	-- 			(self `on` buttonReleaseEvent, onmouseup)] $ \(signal, prop) -> perhaps (prop eventmap) (\listener -> signal $ listener stateref)
	where
		simpleBind = maybeBind self eventmap stateref

-- |
pass :: Monad m => m ()
pass = return ()

-- |
-- bind :: WidgetClass self => self -> Signal ->
bind :: WidgetClass self => self -> Signal self callback -> IORef s -> (IORef s -> callback) -> IO (ConnectId self)
bind self event stateref action = self `on` event $ action stateref

-- |
maybeBind :: WidgetClass self => self -> EventMap self s -> IORef s -> Signal self callback -> (EventMap self s -> Maybe (IORef s -> callback)) -> IO ()
maybeBind self eventmap state event find = maybe pass (void . bind self event state) (find eventmap)

-- |
makeWidget :: WidgetClass self => IO self -> IORef s -> WidgetSettings self s -> IO self
makeWidget make stateref settings = do
	widget <- make
	applySettings widget stateref settings

-- |
applySettings :: WidgetClass self => self -> IORef s -> WidgetSettings self s -> IO self
applySettings self stateref settings = do
	widgetSetSizeRequest self dx dy -- TODO: Use windowSetDefaultSize instead (?)
	widgetAddEvents self eventmasks
	attachListeners self stateref (_settinglisteners settings)
	return self
	where
		(dx:+dy)   = _settingsize settings
		eventmasks = _settingeventmasks settings

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

-- |
-- TODO: Add setup argument (?)
-- TODO: Use 'IO (App s)' or 'App s' directly (?)
run :: IO (App layout state) -> IO ()
run makeApp = makeApp >> mainGUI

max_stars_repo_name: SwiftsNamesake/Southpaw

max_issues_repo_name: SwiftsNamesake/Southpaw

max_forks_repo_name: SwiftsNamesake/Southpaw


text:
module Nondet where
import Data.Complex
import Test.QuickCheck

type Nondet a = [a]

sqrts :: Floating a => a -> Nondet a
sqrts x = [y, negate y]
  where
    y = sqrt x

-- Exercise Four:  Apply the function on all input values and aggregate the results
bind :: (a -> Nondet b) -> (Nondet a -> Nondet b)
bind = undefined

-- Using # instead of * for composition to avoid ambiguities
(#) :: (b -> Nondet c) -> (a -> Nondet b) -> (a -> Nondet c)
g' # f' = bind g' . f'

-- Exercise Five: Provide minimal context for the given value
unit :: a -> Nondet a
unit = undefined

-- lift --- lifting functions
lift :: (a -> b) -> (a -> Nondet b)
lift f = unit . f

-- Solution to the quadratic equation a*x^2 + b * y + c = 0
-- delta = b^2 - 4*a*c
-- x = (-b +- sqrt delta) / (2*a)
solveQuadratic :: Floating a  => a -> a -> a -> Nondet a
solveQuadratic a b c = lift (/(2*a)) # lift (b+) # sqrts $ delta
  where
    delta = b*b - 4*a*c

-- Exercise Six:  Test that (for a given value x)
-- (a) f # unit = unit # f = f
-- (b) lift g # lift f = lift (g.f)
check_unit1, check_unit2 :: (Complex Float -> Nondet (Complex Float)) -> Complex Float -> Bool
check_unit1 f x = undefined
check_unit2 f x = undefined

test_unit1, test_unit2 :: IO ()
test_unit1 = quickCheck $ check_unit1 sqrts
test_unit2 = quickCheck $ check_unit2 sqrts

check_lift :: (Float -> Float) -> (Float -> Float) -> Float -> Bool
check_lift f g x = undefined

test_lift :: IO ()
test_lift = quickCheck $ check_lift (+2) (*3)

-- Exercise Ten(b): Rewrite the module to make Nondet instance of
-- the Monad typeclass
-- Note: You first need to make it an instance of Functor and Applicative

-- Exercise Twelve:  Write the solution to the quadratic equation in do notation

max_stars_repo_name: PavelClaudiuStefan/FMI

max_issues_repo_name: PavelClaudiuStefan/FMI

max_forks_repo_name: PavelClaudiuStefan/FMI


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

module Neural(
             Example,
             initNetwork,
             train,
             predict
             ) where

  import Numeric.LinearAlgebra
  import System.Random
  import Debug.Trace
  import Data.List

  import Tool

  type Example = (Vector Double, Vector Double)

  -- Fully connected neural network described in ADT
  data Network where
    -- Start layer : data layer. consists of just one data vector. (homogeneous)
    StartLayer :: Vector Double -> Network
    -- Hidden layer : consists of three elements
    -- 1. network subset that comes before this layer
    -- 2. weight matrix of size ((curr hidden layer dimension), (prev layer output dimension+1)).
    -- 3. current forward values (homogeneous)
    HiddenLayer :: Network -> Matrix Double -> Vector Double -> Network
    -- End layer : softmax layer. consists of three elements
    -- 1. network subset that comes before this layer
    -- 2. current softmax value
    -- 3. current loss of training (not used)
    EndLayer :: Network -> Vector Double -> Double -> Network
    deriving (Show)

  {- Activations and loss -}

  softmax :: Vector Double -> Vector Double
  softmax v = normalizer `seq` ((\val -> val / normalizer) `cmap` (exp v))
    where normalizer = sumElements (exp v)

  crossEntropy :: Vector Double -> Vector Double -> Double
  crossEntropy t y = -(sumElements $ (t * (log y)))

  activation x = tanh x

  {- Derivation computations -}

  {-
   - jaNodeHidden is a matrix of derivatives
   -
   - dy1/dx1 dy1/dx2 ... dy1/dxr
   -          ...
   - dyk/dx1 dyk/dx2 ... dyk/dxr
   -
   - where yi = activation (sum (weight * xj)), dimension of x is r, dimension of y is k.
   -}
  jaNodeHidden :: Matrix Double -> Vector Double -> Matrix Double
  jaNodeHidden weights y = weights * (repeatCols (cols weights) jaActivation)
    where
      jaActivation = 1 - (y^2)

  {-
   - jaWeightHidden is a matrix of derivatives
   -
   - dy1/dw11 dy1/dw12 ... dy1/dw1r
   -          ...
   - dyk/dwk1 dyk/dwk2 ... dyk/dwkr
   -
   - where yi = activation (sum (wij * xj)), dimension of x is r, dimension of y is k.
   -}
  jaWeightHidden :: Vector Double -> Vector Double -> Matrix Double
  jaWeightHidden x y = jaActivation `outer` jaLinear
    where
      jaActivation = 1 - (y^2)
      jaLinear = x

  -- Derivatives of softmax-crossentropy loss function
  jaLoss t y = y - t

  {- Network handling -}

  -- Initialize network with random weights and zero values
  -- The first argument is dimension of hidden layers.
  --   0th element represents dimension of the last hidden layer.
  --   0th element should be equal to the number of classes
  -- The second argument is dimension of input layer.
  -- length of hiddenDims should be larger than or equal to 1
  initNetwork :: [Int] -> Int -> IO Network
  initNetwork hiddenDims inputDim = do
    inner <- initInnerNetwork hiddenDims inputDim
    return $ EndLayer
      inner
      (zeros $ head hiddenDims)
      0.0
      where
        initInnerNetwork :: [Int] -> Int -> IO Network
        initInnerNetwork hiddenDims inputDim =
          if 1 == length hiddenDims then do
            randWeight <- randmat (-0.5, 0.5) (head hiddenDims) (inputDim+1) -- col+1 for bias
            return $ HiddenLayer
              (StartLayer (homo $ ones inputDim))
              randWeight
              (homo $ zeros (head hiddenDims)) -- append 1 for bias
          else do
            randWeight <- randmat (-0.5, 0.5) (head hiddenDims) (head (tail hiddenDims) + 1) -- col+1 for bias
            prevLayer <- (initInnerNetwork (tail hiddenDims) inputDim)
            return $ HiddenLayer
              prevLayer
              randWeight
              (homo $ zeros (last hiddenDims)) -- append 1 for bias

  -- Get a vector value of the last layer of the network
  lastLayerOf :: Network -> Vector Double
  lastLayerOf network = case network of
    StartLayer d -> d
    HiddenLayer _ _ value -> value
    EndLayer _ value _ -> value

  -- Train FC neural network 'epoch' times using given examples
  train :: Int -> Double -> [Example] -> Network -> Network
  train epoch learningRate examples network = examples `seq` ((iterate (trainEpoch learningRate examples) network)!!epoch)

  -- Train FC neural network with examples
  trainEpoch :: Double -> [Example] -> Network -> Network
  trainEpoch learningRate examples network =
    foldl' (\acc (input, gt) ->
      ((backward learningRate gt) . forward . (feedData input)) acc
      ) network examples

  -- Predict softmax values of network
  predict :: Network -> Vector Double -> Vector Double
  predict network input = (lastLayerOf . forward . (feedData input)) network

  -- Prepare Network with input
  feedData :: Vector Double -> Network -> Network
  feedData input network = input `seq` case network of
    StartLayer d -> StartLayer $ homo input
    HiddenLayer prevLayer weights value -> HiddenLayer (feedData input prevLayer) weights value
    EndLayer prevLayer value loss -> EndLayer (feedData input prevLayer) value loss

  -- Forward data to update values of Network. Leaves weights of layers unchanged.
  -- Use this for prediction
  forward :: Network -> Network
  forward network = case network of
    StartLayer d -> StartLayer d
    HiddenLayer prevLayer weights _ ->
      value `seq` HiddenLayer prevResult weights value
        where prevResult = forward prevLayer
              value = (homo (activation (weights #> lastLayerOf prevResult)))
    EndLayer prevLayer _ loss ->
      lastVector `seq` EndLayer prevResult (softmax $ invhomo lastVector) loss -- remove the last element since it is meaningless
        where prevResult = forward prevLayer
              lastVector = lastLayerOf prevResult

  -- Backpropagate errors to update weights of Network. Leaves values of layers unchanged
  -- Uses stochastic gradient descent
  backward :: Double -> Vector Double -> Network -> Network
  backward learningRate gt (EndLayer prevLayer value _) =
    EndLayer (innerBackward prevLayer lastDiff) value loss
      where
        {-
         - lastDiff is a vector of partial derivatives
         -
         - [dL/dY1 ... dL/dYk]
         -
         - where last softmax layer is k-dimension.
         -}
        lastDiff = jaLoss gt value
        loss = crossEntropy gt value

        innerBackward :: Network -> Vector Double -> Network
        innerBackward network nextDiffs = case network of
          StartLayer d -> StartLayer d
          HiddenLayer prevLayer weights values -> newWeights `seq` HiddenLayer (innerBackward prevLayer currDiff) newWeights values
            where
              {-
               - currDiff is a vector of partial derivatives
               -
               - [dL/dX1 ... dL/dXr]
               -
               - where X is the output of previous layer
               -}
              currDiff = (invhomo . flatten) $ ((jaNodeHidden weights (invhomo values)) * (repeatCols (cols weights) nextDiffs)) ? [0]
              newWeights = weightDiff `seq` (weights - ((\w -> w*learningRate) `cmap` weightDiff))
                where
                  {-
                   - weightDiff is a matrix filled with partial derivatives
                   -
                   -  dL/dW11  dL/dW12 ...  dL/dW1r
                   -               ...
                   -  dL/dWk1  dL/dWk2 ...  dL/dWkr
                   -
                   - where this hidden layer have 'k' nodes and previous layer has 'r' nodes
                   -}
                  weightDiff = pjacobian `seq` (pjacobian * (repeatCols (cols pjacobian) nextDiffs))
                    where
                      {-
                       - pjacobian is a matrix filled with partial derivatives (p is for pseudo)
                       -
                       -  dY1/dW11  dY1/dW12 ...  dY1/dW1r
                       -               ...
                       -  dYk/dWk1  dYk/dWk2 ...  dYk/dWkr
                       -
                       - where this hidden layer have 'k' nodes and previous layer has 'r' nodes
                       -}
                      pjacobian = let x = lastLayerOf prevLayer in x `seq` jaWeightHidden x (invhomo values)

max_stars_repo_name: ycjungSubhuman/hskmlp

max_issues_repo_name: ycjungSubhuman/hskmlp

max_forks_repo_name: ycjungSubhuman/hskmlp


text:
module NeuralNet where

import MatrixPlus

import Numeric.LinearAlgebra as N
import Test.QuickCheck

class NeuralNet n where
  evaluate
    :: n        -- ^ The neural net
    -> [Double] -- ^ The input pattern
    -> [Double] -- ^ The output pattern
  train
    :: n        -- ^ The neural net before training
    -> [Double] -- ^ The input pattern
    -> [Double] -- ^ The target pattern
    -> n        -- ^ The neural net after training

-- {{{ Useful things for testing

-- | Generate an arbitrary number between zero and one.
-- | To see sample values, in GHCi type: sample arbZeroToOne
arbZeroToOne :: Gen Double
arbZeroToOne = choose (0, 1)

-- | Generate a column vector of the specified length with arbitrary values.
-- | To see sample values (of length 4), in GHCi type: sample (sizedArbColumnVector 4)
sizedArbColumnVector :: Int -> Gen (ColumnVector Double)
sizedArbColumnVector n = do
    values <- vectorOf n arbZeroToOne
    return (listToColumnVector values)

-- | Generate a column vector where all values are zero
zeroColumnVector :: Int -> ColumnVector Double
zeroColumnVector n = listToColumnVector (replicate n 0)

-- | Generate a weight matrix of the specified dimensions with arbitrary values.
-- | To see sample values (of size 3x2), in GHCi type: sample (sizedArbWeightMatrix 3 2)
sizedArbWeightMatrix :: Int -> Int -> Gen (Matrix Double)
sizedArbWeightMatrix r c = do
    values <- vectorOf (r*c) arbZeroToOne
    return ((r N.><c) values)

-- }}}

max_stars_repo_name: mhwombat/backprop-example

max_issues_repo_name: mhwombat/backprop-example

max_forks_repo_name: mhwombat/backprop-example


text:
{-# OPTIONS -cpp  #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Control.Parallel.Strategies
-- Copyright   :  (c) The University of Glasgow 2001
-- License     :  BSD-style (see the file libraries/base/LICENSE)
--
-- Maintainer  :  libraries@haskell.org
-- Stability   :  experimental
-- Portability :  non-portable
--
-- Parallel strategy combinators. See
-- <http://www.macs.hw.ac.uk/~dsg/gph/papers/html/Strategies/strategies.html>
-- for more information.
--
-- Original authors:
--	Phil Trinder, Hans-Wolfgang Loidl, Kevin Hammond et al.
--
-----------------------------------------------------------------------------
-- #ifdef __PARALLEL_HASKELL__
-- module Control.Parallel.Strategies (
-- #else
module Parallel.Strategies (
-- #endif
   -- *	Strategy Type, Application and Semantics
   Done, Strategy,
   (>|), (>||),
   using, demanding, sparking,
   -- *	Basic Strategies
   r0, rwhnf, NFData(..),
   -- * Strategic Function Application
   ($|), ($||),
   (.|), (.||),
   (-|), (-||),
   -- * Tuples
   seqPair, parPair,
   seqTriple, parTriple,
   -- * Lists: Parallel Strategies
   parList, parListN, parListNth, parListChunk,
   parMap, parFlatMap, parZipWith,
   -- * Lists: Sequential Strategies
   seqList, seqListN, seqListNth, parBuffer,
   -- *	Arrays
   seqArr, parArr,
   -- * Deprecated types and functions
   sPar, sSeq,
   Assoc(..),
   fstPairFstList, force, sforce
  ) where

-- based on hslibs/concurrent/Strategies.lhs; see it for more detailed
-- code comments.
#ifdef __PARALLEL_HASKELL__
#else
#warning "Using nonparallel GHC. The following line will fail."
#endif
import Control.Parallel as Parallel (pseq)
import Data.Array
import Data.Complex
import Data.Int
import qualified Data.IntMap (IntMap, toList)
import qualified Data.IntSet (IntSet, toList)
import qualified Data.Map (Map, toList)
import qualified Data.Set (Set, toList)
import qualified Data.Tree (Tree(..))
import Data.Word

import Prelude hiding (seq)
import qualified Prelude (seq)

-- not a terribly portable way of getting at Ratio rep.
#ifdef __GLASGOW_HASKELL__
import GHC.Real	(Ratio(..))	-- The basic defns for Ratio
#endif

#ifdef __HUGS__
import Hugs.Prelude(Ratio(..) )
#endif

#ifdef __NHC__
import Ratio (Ratio(..) )
#endif

infixl 0 `using`,`demanding`,`sparking`              -- weakest precedence!

infixr 2 >||                -- another name for par
infixr 3 >|                 -- another name for seq
infixl 6 $||, $|            -- strategic function application (seq and par)
infixl 9 .|, .||, -|, -||   -- strategic (inverse) function composition

-- We need 'pseq', not the Prelude 'seq' here.  See the documentation
-- with 'pseq' in Control.Parallel.
seq = Parallel.pseq
par x y = y

------------------------------------------------------------------------------
-- *			Strategy Type, Application and Semantics
------------------------------------------------------------------------------

{-
The basic combinators for strategies are 'par' and 'seq' but with types that
indicate that they only combine the results of a strategy application.

NB: This version can be used with Haskell 1.4 (GHC 2.05 and beyond), *but*
    you won't get strategy checking on seq (only on par)!

The operators >| and >|| are alternative names for `seq` and `par`.
With the introduction of a Prelude function `seq` separating the Prelude
function from the Strategy function becomes a pain. The notation also matches
the notation for strategic function application.
-}

type Done = ()

-- | A strategy takes a value and returns a 'Done' value to indicate that
--   the specifed evaluation has been performed.
type Strategy a = a -> Done

-- | Evaluates the first argument before the second.
(>|) :: Done -> Done -> Done
{-# INLINE (>|) #-}
(>|) = Prelude.seq

-- | Evaluates the first argument in parallel with the second.
(>||) :: Done -> Done -> Done
{-# INLINE (>||) #-}
(>||) = par -- Parallel.par

-- | Takes a value and a strategy, and applies the strategy to the
-- value before returning the value. Used to express data-oriented
-- parallelism. @x \`using\` s@ is a projection on @x@, i.e. both:
--
-- [a retraction] @x \`using\` s@ ⊑ @x@
--
-- [idempotent] @(x \`using\` s) \`using\` s@ = @x \`using\` s@
--
using :: a -> Strategy a -> a
using x s = s x `seq` x

-- | Evaluates the second argument before the first.
-- Used to express control-oriented parallelism. The second
-- argument is usually a strategy application.
demanding :: a -> Done -> a
demanding = flip seq

-- | Evaluates the second argument in parallel with the first.
-- Used to express control-oriented
-- parallelism. The second argument is usually a strategy application.
sparking :: a -> Done -> a
sparking  = flip par -- Parallel.par
-- Sparking should only be used
-- with a singleton sequence as it is not necessarily executed.

-- | A strategy corresponding to 'par':
-- @x \`par\` e@ = @e \`using\` sPar x@.
--
-- 'sPar' has been superceded by 'sparking'.
-- Replace @e \`using\` sPar x@ with @e \`sparking\` rwhnf x@.
{-# DEPRECATED sPar "Use sparking instead." #-}
sPar :: a -> Strategy b
sPar x y = x `par` ()

-- | A strategy corresponding to 'seq':
-- @x \`seq\` e@ = @e \`using\` sSeq x@.
--
-- 'sSeq' has been superceded by 'demanding'.
-- Replace @e \`using\` sSeq x@ with @e \`demanding\` rwhnf x@.
{-# DEPRECATED sSeq "Use demanding instead." #-}
sSeq :: a -> Strategy b
sSeq x y = x `seq` ()

-----------------------------------------------------------------------------
-- *			Basic Strategies
-----------------------------------------------------------------------------

-- | Performs /no/ evaluation of its argument.
r0 :: Strategy a
r0 x = ()

-- | Reduces its argument to weak head normal form.
rwhnf :: Strategy a
rwhnf x = x `seq` ()

class NFData a where
  -- | Reduces its argument to (head) normal form.
  rnf :: Strategy a
  -- Default method. Useful for base types. A specific method is necessay for
  -- constructed types
  rnf = rwhnf

class (NFData a, Integral a) => NFDataIntegral a
class (NFData a, Ord a) => NFDataOrd a

------------------------------------------------------------------------------
-- *                     Strategic Function Application
------------------------------------------------------------------------------

{-
These are very
handy when writing pipeline parallelism asa sequence of @$@, @$|@ and
@$||@'s. There is no need of naming intermediate values in this case. The
separation of algorithm from strategy is achieved by allowing strategies
only as second arguments to @$|@ and @$||@.
-}

-- | Sequential function application. The argument is evaluated using
--   the given strategy before it is given to the function.
($|) :: (a -> b) -> Strategy a -> a -> b
f $| s  = \ x -> f x `demanding` s x

-- | Parallel function application. The argument is evaluated using
-- the given strategy, in parallel with the function application.
($||) :: (a -> b) -> Strategy a -> a -> b
f $|| s = \ x -> f x `sparking` s x

-- | Sequential function composition. The result of
-- the second function is evaluated using the given strategy,
-- and then given to the first function.
(.|) :: (b -> c) -> Strategy b -> (a -> b) -> (a -> c)
(.|) f s g = \ x -> let  gx = g x
                    in   f gx `demanding` s gx

-- | Parallel function composition. The result of the second
-- function is evaluated using the given strategy,
-- in parallel with the application of the first function.
(.||) :: (b -> c) -> Strategy b -> (a -> b) -> (a -> c)
(.||) f s g = \ x -> let  gx = g x
                     in   f gx `sparking` s gx

-- | Sequential inverse function composition,
-- for those who read their programs from left to right.
-- The result of the first function is evaluated using the
-- given strategy, and then given to the second function.
(-|) :: (a -> b) -> Strategy b -> (b -> c) -> (a -> c)
(-|) f s g = \ x -> let  fx = f x
                    in   g fx `demanding` s fx

-- | Parallel inverse function composition,
-- for those who read their programs from left to right.
-- The result of the first function is evaluated using the
-- given strategy, in parallel with the application of the
-- second function.
(-||) :: (a -> b) -> Strategy b -> (b -> c) -> (a -> c)
(-||) f s g = \ x -> let  fx = f x
                     in   g fx `sparking` s fx

------------------------------------------------------------------------------
--			Marking a Strategy
------------------------------------------------------------------------------

{-
Marking a strategy.

Actually, @markStrat@  sticks a label @n@  into the sparkname  field of the
thread executing strategy @s@. Together with a runtime-system that supports
propagation of sparknames to the children this means that this strategy and
all its children have  the sparkname @n@ (if the  static sparkname field in
the @parGlobal@ annotation contains the value 1). Note, that the @SN@ field
of starting the marked strategy itself contains the sparkname of the parent
thread. The END event contains @n@ as sparkname.
-}

#if 0
markStrat :: Int -> Strategy a -> Strategy a
markStrat n s x = unsafePerformPrimIO (
     _casm_ ``%r = set_sparkname(CurrentTSO, %0);'' n `thenPrimIO` \ z ->
     returnPrimIO (s x))
#endif

-----------------------------------------------------------------------------
--			Strategy Instances and Functions
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
-- *	                Tuples
-----------------------------------------------------------------------------

{-
We currently support up to 9-tuples. If you need longer tuples you have to
add the instance explicitly to your program.
-}

instance (NFData a, NFData b) => NFData (a,b) where
  rnf (x,y) = rnf x `seq` rnf y

instance (NFData a, NFData b, NFData c) => NFData (a,b,c) where
  rnf (x,y,z) = rnf x `seq` rnf y `seq` rnf z

instance (NFData a, NFData b, NFData c, NFData d) => NFData (a,b,c,d) where
  rnf (x1,x2,x3,x4) = rnf x1 `seq`
		        rnf x2 `seq`
		        rnf x3 `seq`
		        rnf x4

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5) =>
         NFData (a1, a2, a3, a4, a5) where
  rnf (x1, x2, x3, x4, x5) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6) =>
         NFData (a1, a2, a3, a4, a5, a6) where
  rnf (x1, x2, x3, x4, x5, x6) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6, NFData a7) =>
         NFData (a1, a2, a3, a4, a5, a6, a7) where
  rnf (x1, x2, x3, x4, x5, x6, x7) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6 `seq`
                  rnf x7

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6, NFData a7, NFData a8) =>
         NFData (a1, a2, a3, a4, a5, a6, a7, a8) where
  rnf (x1, x2, x3, x4, x5, x6, x7, x8) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6 `seq`
                  rnf x7 `seq`
                  rnf x8

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6, NFData a7, NFData a8, NFData a9) =>
         NFData (a1, a2, a3, a4, a5, a6, a7, a8, a9) where
  rnf (x1, x2, x3, x4, x5, x6, x7, x8, x9) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6 `seq`
                  rnf x7 `seq`
                  rnf x8 `seq`
                  rnf x9

-- | Apply two strategies to the elements of a pair sequentially
--   from left to right.
seqPair :: Strategy a -> Strategy b -> Strategy (a,b)
seqPair strata stratb (x,y) = strata x `seq` stratb y

-- | Apply two strategies to the elements of a pair in parallel.
parPair :: Strategy a -> Strategy b -> Strategy (a,b)
parPair strata stratb (x,y) = strata x `par` stratb y `par` ()
-- The reason for the last 'par' is so that the strategy terminates
-- quickly. This is important if the strategy is used as the 1st
-- argument of a seq

-- | Apply three strategies to the elements of a triple in sequentially
--   from left to right.
seqTriple :: Strategy a -> Strategy b -> Strategy c -> Strategy (a,b,c)
seqTriple strata stratb stratc p@(x,y,z) =
  strata x `seq`
  stratb y `seq`
  stratc z

-- | Apply three strategies to the elements of a triple in parallel.
parTriple :: Strategy a -> Strategy b -> Strategy c -> Strategy (a,b,c)
parTriple strata stratb stratc (x,y,z) =
  strata x `par`
  stratb y `par`
  stratc z `par`
  ()

-----------------------------------------------------------------------------
-- 			Atomic types
-----------------------------------------------------------------------------

{-
Weak head normal form and normal form are identical for integers, so the
default rnf is sufficient.
-}
instance NFData Int
instance NFData Integer
instance NFData Float
instance NFData Double

instance NFData Int8
instance NFData Int16
instance NFData Int32
instance NFData Int64

instance NFData Word8
instance NFData Word16
instance NFData Word32
instance NFData Word64

instance NFDataIntegral Int
instance NFDataOrd Int

--Rational and complex numbers.

instance (Integral a, NFData a) => NFData (Ratio a) where
  rnf (x:%y) = rnf x `seq`
               rnf y `seq`
               ()

instance (RealFloat a, NFData a) => NFData (Complex a) where
  rnf (x:+y) = rnf x `seq`
	         rnf y `seq`
               ()

instance NFData Char
instance NFData Bool
instance NFData ()

-----------------------------------------------------------------------------
-- 			Various library types
-----------------------------------------------------------------------------

instance NFData a => NFData (Maybe a) where
    rnf Nothing  = ()
    rnf (Just x) = rnf x

instance (NFData a, NFData b) => NFData (Either a b) where
    rnf (Left x)  = rnf x
    rnf (Right y) = rnf y

instance (NFData k, NFData a) => NFData (Data.Map.Map k a) where
    rnf = rnf . Data.Map.toList

instance NFData a => NFData (Data.Set.Set a) where
    rnf = rnf . Data.Set.toList

instance NFData a => NFData (Data.Tree.Tree a) where
    rnf (Data.Tree.Node r f) = rnf r `seq` rnf f

instance NFData a => NFData (Data.IntMap.IntMap a) where
    rnf = rnf . Data.IntMap.toList

instance NFData Data.IntSet.IntSet where
    rnf = rnf . Data.IntSet.toList

-----------------------------------------------------------------------------
-- 			Lists
-----------------------------------------------------------------------------

instance NFData a => NFData [a] where
  rnf [] = ()
  rnf (x:xs) = rnf x `seq` rnf xs

----------------------------------------------------------------------------
-- *                   Lists: Parallel Strategies
----------------------------------------------------------------------------

-- | Applies a strategy to every element of a list in parallel.
parList :: Strategy a -> Strategy [a]
parList strat []     = ()
parList strat (x:xs) = strat x `par` (parList strat xs)

-- | Applies a strategy to the first @n@ elements of a list in parallel.
parListN :: (Integral b) => b -> Strategy a -> Strategy [a]
parListN n strat []     = ()
parListN 0 strat xs     = ()
parListN n strat (x:xs) = strat x `par` (parListN (n-1) strat xs)

-- | Evaluates @n@ elements of the spine of the argument list and applies
-- the given strategy to the @n@th element (if there is one) in parallel with
-- the result. E.g. @parListNth 2 [e1, e2, e3]@ evaluates @e3@.
parListNth :: Int -> Strategy a -> Strategy [a]
parListNth n strat xs
  | null rest = ()
  | otherwise = strat (head rest) `par` ()
  where
    rest = drop n xs

-- | Splits a list into chunks (sub-sequences) of length @n@,
-- and applies a strategy sequentially to the elements in each
-- chunk. The chunks are evaluated in parallel.
-- This is useful for increasing the grain size.
parListChunk :: Int -> Strategy a -> Strategy [a]
parListChunk n strat [] = ()
parListChunk n strat xs = seqListN n strat xs `par`
			    parListChunk n strat (drop n xs)

-- | Applies a function to each element of a list and
-- and evaluates the result list in parallel,
-- using the given strategy for each element.
parMap :: Strategy b -> (a -> b) -> [a] -> [b]
parMap strat f xs = map f xs `using` parList strat

-- | Uses 'parMap' to apply a list-valued function to each
-- element of a list in parallel, and concatenates the results.
parFlatMap :: Strategy [b] -> (a -> [b]) -> [a] -> [b]
parFlatMap strat f xs = concat (parMap strat f xs)

-- | Zips together two lists using a function,
-- and evaluates the result list in parallel.
parZipWith :: Strategy c -> (a -> b -> c) -> [a] -> [b] -> [c]
parZipWith strat z as bs =
  zipWith z as bs `using` parList strat

----------------------------------------------------------------------------
-- *                     Lists: Sequential Strategies
----------------------------------------------------------------------------

-- | Sequentially applies a strategy to each element of a list.
seqList :: Strategy a -> Strategy [a]
seqList strat []     = ()
seqList strat (x:xs) = strat x `seq` (seqList strat xs)

-- | Sequentially applies a strategy to the first n elements of a list.
seqListN :: (Integral a) => a -> Strategy b -> Strategy [b]
seqListN n strat []     = ()
seqListN 0 strat xs     = ()
seqListN n strat (x:xs) = strat x `seq` (seqListN (n-1) strat xs)

-- | Applies a strategy to the @n@th element of a list
--  (if there is one) before returning the result.
--  E.g. @seqListNth 2 [e1, e2, e3]@ evaluates @e3@.
seqListNth :: Int -> Strategy b -> Strategy [b]
seqListNth n strat xs
  | null rest = ()
  | otherwise = strat (head rest)
  where
    rest = drop n xs

-- | Parallel n-buffer function added for the revised version of the strategies
-- paper. 'parBuffer' supersedes the older @fringeList@. It has the same
-- semantics.
parBuffer :: Int -> Strategy a -> [a] -> [a]
parBuffer n s xs =
  return xs (start n xs)
  where
    return (x:xs) (y:ys) = (x:return xs ys) `sparking` s y
    return xs     []     = xs

    start n []     = []
    start 0 ys     = ys
    start n (y:ys) = start (n-1) ys `sparking` s y

{-
 'fringeList' implements a `rolling buffer' of length n, i.e.applies a
 strategy to the nth element of list when the head is demanded. More
 precisely:

   semantics:         fringeList n s = id :: [b] -> [b]
   dynamic behaviour: evalutates the nth element of the list when the
		      head is demanded.

 The idea is to provide a `rolling buffer' of length n.
fringeList :: (Integral a) => a -> Strategy b -> [b] -> [b]
fringeList n strat [] = []
fringeList n strat (r:rs) =
  seqListNth n strat rs `par`
  r:fringeList n strat rs
-}

------------------------------------------------------------------------------
-- *			Arrays
------------------------------------------------------------------------------
instance (Ix a, NFData a, NFData b) => NFData (Array a b) where
  rnf x = rnf (bounds x) `seq` seqList rnf (elems x) `seq` ()

-- | Apply a strategy to all elements of an array sequentially.
seqArr :: (Ix b) => Strategy a -> Strategy (Array b a)
seqArr s arr = seqList s (elems arr)

-- | Apply a strategy to all elements of an array in parallel.
parArr :: (Ix b) => Strategy a -> Strategy (Array b a)
parArr s arr = parList s (elems arr)

{-# DEPRECATED Assoc "Does not belong in Control.Parallel.Strategies" #-}
data  Assoc a b =  a := b  deriving ()

instance (NFData a, NFData b) => NFData (Assoc a b) where
  rnf (x := y) = rnf x `seq` rnf y `seq` ()

------------------------------------------------------------------------------
-- *	                Some strategies specific for Lolita
------------------------------------------------------------------------------

{-# DEPRECATED fstPairFstList "This was just an example. Write your own." #-}
fstPairFstList :: (NFData a) => Strategy [(a,b)]
fstPairFstList = seqListN 1 (seqPair rwhnf r0)

-- Some HACKs for Lolita. AFAIK force is just another name for our rnf and
-- sforce is a shortcut (definition here is identical to the one in Force.lhs)

{-# DEPRECATED force, sforce "Lolita-specific hacks." #-}
force :: (NFData a) => a -> a
sforce :: (NFData a) => a -> b -> b

force = id $| rnf
sforce x y = force x `seq` y

max_stars_repo_name: Parrows/Parrows

max_issues_repo_name: Parrows/Parrows

max_forks_repo_name: Parrows/Parrows


text:
{-# LANGUAGE RankNTypes #-}

module Kalman
  ( movingAverage,
    System (..),
    Estimate (..),
    Measurement (..),
    genSigmas,
    ekf,
    ut,
    ukf,
    weightedCrossCov,
    runKF,
  )
where

import Control.Applicative
import Numeric.AD
import Numeric.AD.Internal.Reverse
import Numeric.LinearAlgebra
import Sim
import System.Random

type Sigmas = [[Double]]

data System = System
  { fA, fH :: forall a. Floating a => ([a] -> [a]),
    kQ, kR :: Matrix Double
  }

data Estimate = Estimate
  { sX :: Vector Double,
    sP :: Matrix Double
  }
  deriving (Show)

-- Simple 2 tap FIR to appx moving average filter
movingAverage :: [Vector Double] -> [Vector Double]
movingAverage = scanl1 (\v1 v2 -> (0.9 * v1) + (0.1 * v2))

-- kalman :: System -> Estimate -> Measurement -> Estimate
-- kalman (System a c q r) (Estimate x p) z = Estimate x' p' where
--   -- Convert lists to matrices since we have linear system
--   lA = fromLists a
--   lC = fromLists c
--   -- prediction step
--   px = lA #> x
--   pp = lA <> p <> tr lA + q
--   -- update step
--   g = pp <> tr lC <> inv (lC <> pp <> tr lC + r)
--   x' = px + g #> (z - lC #> px)
--   p' = (ident (size x) - g <> lC) <> pp

runKF :: (System -> Estimate -> Measurement -> Estimate) -> System -> Estimate -> [Measurement] -> [Estimate]
runKF kf sys = scanl (kf sys)

ekf :: System -> Estimate -> Measurement -> Estimate
ekf (System fA fH kQ kR) (Estimate x p) z = Estimate x' p'
  where
    -- prediction step
    a = fromLists (jacobian fA (toList x))
    px = fA $ toList x
    pp = a <> p <> tr a + kQ
    -- update step
    h = fromLists $ jacobian fH px
    k = pp <> tr h <> inv (h <> pp <> tr h + kR)
    x' = vector px + k #> (z - vector (take (size z) $ fH px))
    p' = (ident (size x) - k <> h) <> pp

genSigmas :: Estimate -> Double -> Sigmas
genSigmas (Estimate x p) c = sigmas
  where
    -- Cholesky Factorization to get sqrt of matrix
    cP = scale (sqrt c) $ chol $ sym p
    -- Add/sub the sqrt covariance matrix from the current estimate of x
    xs = fromLists $ replicate (size x) (toList x)
    pos = toLists $ xs + cP
    neg = toLists $ xs - cP
    -- Combine all sigma points including the original estimate of x at the head
    sigmas = [toList x] ++ pos ++ neg

-- Unscented Transformation
ut :: [Vector Double] -> [Double] -> [Double] -> Estimate
ut xs wm wc = Estimate x' p'
  where
    -- Weight each sigma by scaling factor
    xw = getZipList $ scale <$> ZipList wm <*> ZipList xs
    x' = sum xw
    -- Compute forecasted state statistics
    v1 = map (\v -> fromColumns [v - x']) xs
    v2 = map (\v -> fromRows [v - x']) xs
    cov = getZipList $ (<>) <$> ZipList v1 <*> ZipList v2
    covw = getZipList $ scale <$> ZipList wc <*> ZipList cov
    p' = sum covw

weightedCrossCov :: Vector Double -> [Vector Double] -> Vector Double -> [Vector Double] -> [Double] -> Matrix Double
weightedCrossCov v1 vs1 v2 vs2 wc = cP
  where
    vs1' = map (\v -> fromColumns [v - v1]) vs1
    vs2' = map (\v -> fromRows [v - v2]) vs2
    cov = getZipList $ (<>) <$> ZipList vs1' <*> ZipList vs2'
    covw = getZipList $ scale <$> ZipList wc <*> ZipList cov
    cP = sum covw

ukf :: System -> Estimate -> Measurement -> Estimate
ukf (System fA fH kQ kR) (Estimate x p) z = Estimate xh ph
  where
    l = fromIntegral $ size x
    m = fromIntegral $ size z
    alpha = 1e-3
    ki = 0
    beta = 2

    lambda = alpha ^ 2 * (l + ki) - l
    c = l + lambda

    wm = (lambda / c) : replicate (2 * fromIntegral (size x)) (0.5 / c)
    wc = ((lambda / c) + (1 - alpha ^ 2 + beta)) : replicate (2 * fromIntegral (size x)) (0.5 / c)

    -- Generate sigma points around current estimate
    sigmas = genSigmas (Estimate x p) c

    -- Run sigmas through non linear function to produce new points
    xs = map (fromList . fA) sigmas
    zs = map (fromList . fH) sigmas

    Estimate x' cX = ut xs wm wc
    Estimate z' cZ = ut zs wm wc

    -- Calculate the weighted cross covariance
    cXZ = weightedCrossCov x' xs z' zs wc

    kK = cXZ <> inv (cZ + kR)
    xh = x' + kK #> (z - z')
    ph = (cX + kQ) - kK <> (cZ + kR) <> tr kK


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

-- Created _ _ 2015

-- TODO | -
--        -

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC Pragmas
--------------------------------------------------------------------------------------------------------------------------------------------

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module Cartesian.Types (
  -- ^ Third party types
  V1(..), V2(..), V3(..), V4(..), Complex(..),

  -- ^ Synonyms
  BoxLens, Axis, Axes, Polygon,

  -- ^ Coordinate types
  Normalised, Absolute,

  -- ^ Types defined in this library
  BoundingBox(..), Line, Linear,

  -- ^ Classes
  HasX, HasY, HasZ) where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Linear.V1
import Linear.V2
import Linear.V3
import Linear.V4

import Data.Complex (Complex(..))

import Cartesian.Internal.Types

max_stars_repo_name: jordanemedlock/Cartesian

max_issues_repo_name: jordanemedlock/Cartesian

max_forks_repo_name: jordanemedlock/Cartesian


text:
module Lib where

import qualified Data.Matrix.MatrixMarket as M

import qualified Numeric.LinearAlgebra.Sparse as LA

mmMatrixToSpMatrix (M.RMatrix (m,n) nz str ll) = LA.fromListSM (m,n) ll

max_stars_repo_name: ocramz/sparse-linear-algebra-IO

max_issues_repo_name: ocramz/sparse-linear-algebra-IO

max_forks_repo_name: ocramz/sparse-linear-algebra-IO


text:
module Physics (module Data.Complex, module Physics) where

import Data.Complex

type Vector = Complex Float
j :: Vector
j = 0:+1

normalize :: RealFloat a => Complex a -> Complex a
normalize (0:+0) = 0:+0
normalize a = a / (sqrt $ a * conjugate a)

toxy :: Complex a -> (a, a)
toxy (x :+ y) = (x, y)

toVector :: Float -> Vector
toVector x = x:+0

max_stars_repo_name: jarys/jampazen

max_issues_repo_name: jarys/jampazen

max_forks_repo_name: jarys/jampazen


text:
{-# LANGUAGE FlexibleContexts #-}

------------------------------------------------
-- |
-- Module    :  Numeric.Matrix.Integral.HNFLLL
-- Copyright :  (c) Jun Yoshida 2019
-- License   :  BSD3
--
-- Compute Hermite normal forms in the LLL-based method.
-- The original "pseudo-code" is found in the paper
--   George Havas, Bohdan S. Majewski & Keith R. Matthews (1998) Extended GCD and Hermite Normal Form Algorithms via Lattice Basis Reduction, Experimental Mathematics, 7:2, 125-136, DOI: 10.1080/10586458.1998.10504362
--
------------------------------------------------

module Numeric.Matrix.Integral.HNFLLL (
  hnfLLLST,
  hnfLLL,
  ) where

import Control.Monad (when, forM_)
import Control.Monad.ST
import Control.Monad.Loops (whileM_)

import Data.STRef

import qualified Data.Maybe as M

import Data.Int

import Numeric.LinearAlgebra ((!))
import qualified Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra.Devel

-- | Type class which guarantees the type can be embedded in the real number field
class (Num t,Eq t,Ord t) => Realizable t where
  toReal :: t -> Double

instance Realizable Double where
  toReal = id

instance Realizable Int64 where
  toReal = fromIntegral

-- | Type class to qualify a type to be used in LLL-based algorithm
class (Realizable t, LA.Indexable (LA.Vector t) t, LA.Container LA.Vector t) => LLLQualified t

instance LLLQualified Double

instance LLLQualified Int64

-- | Data Type to carry resources in the computation of LLL-based algorithm
data HNFLLLResource s t = HNFLLL {
  matA :: STMatrix s t,
  matU :: STMatrix s t,
  lambda :: STMatrix s Double,
  focusC1 :: STRef s (Maybe Int),
  focusC2 :: STRef s (Maybe Int) }

-- | The function which is called Reduce2 in the original "pseudo-code."
reduceLLL :: (LLLQualified t) => Int -> Int -> HNFLLLResource s t -> ST s ()
reduceLLL k i resRef = do
  ma <- freezeMatrix (matA resRef)
  case M.listToMaybe (LA.find (/=0) (ma!i)) of
    (Just j) -> do
      writeSTRef (focusC1 resRef) (Just j)
      when ( (ma!i!j) < 0 ) $ do
        rowOper (SCAL (-1) (Row i) AllCols) (matA resRef)
        rowOper (SCAL (-1) (Row i) AllCols) (matU resRef)
        rowOper (SCAL (-1) (Row (j+1)) (ColRange 0 j)) (lambda resRef)
        rowOper (SCAL (-1) (FromRow (j+2)) (Col (j+1))) (lambda resRef)
    Nothing -> writeSTRef (focusC1 resRef) Nothing
  case M.listToMaybe (LA.find (/=0) (ma!k)) of
    (Just j) -> writeSTRef (focusC2 resRef) (Just j)
    Nothing -> writeSTRef (focusC2 resRef) Nothing
  c1 <- readSTRef (focusC1 resRef)
  q <- fromIntegral.floor <$> do
    case c1 of
      (Just j) ->
        (/) <$> (toReal <$> readMatrix (matA resRef) k j) <*> (toReal <$> readMatrix (matA resRef) i j)
      Nothing -> do
        lki <- readMatrix (lambda resRef) (k+1) (i+1)
        di <- readMatrix (lambda resRef) (i+1) (i+1)
        if 2 * abs lki > di
          then return $ fromIntegral (round (lki / di))
          else return 0
  when (q /= 0) $ do
    rowOper (AXPY (negate q) i k AllCols) (matA resRef)
    rowOper (AXPY (negate q) i k AllCols) (matU resRef)
    rowOper (AXPY (toReal $ negate q) (i+1) (k+1) AllCols) (lambda resRef)

-- | The function which is called Swap1 in the original "pseudo-code."
swapLLL :: (Num t, LA.Element t) => Int -> HNFLLLResource s t -> ST s ()
swapLLL k resRef = do
  rowOper (SWAP k (k-1) AllCols) (matA resRef)
  rowOper (SWAP k (k-1) AllCols) (matU resRef)
  rowOper (SWAP k (k+1) (ColRange 0 (k-1))) (lambda resRef)
  ml <- freezeMatrix (lambda resRef)
  forM_ [(k+2)..(LA.cols ml-1)] $ \i -> do
    writeMatrix (lambda resRef) i k $
      ((ml!i!k)*(ml!(k+1)!k) + (ml!i!(k+1))*(ml!(k-1)!(k-1)))/(ml!k!k)
    writeMatrix (lambda resRef) i (k+1) $
      ((ml!i!k)*(ml!(k+1)!(k+1)) + (ml!i!(k+1))*(ml!(k+1)!k))/(ml!k!k)
  writeMatrix (lambda resRef) k k $
    ((ml!(k+1)!k)*(ml!(k+1)!k) + (ml!(k+1)!(k+1))*(ml!(k-1)!(k-1)))/(ml!k!k)

-- | Recognize whether a swap should occur or not.
shouldSwapST :: (Num t, LA.Element t) => Int -> HNFLLLResource s t -> ST s Bool
shouldSwapST k resRef = do
  ml <- freezeMatrix $ lambda resRef
  mc1 <- readSTRef $ focusC1 resRef
  mc2 <- readSTRef $ focusC2 resRef
  case (mc1,mc2) of
    (Just c1,Just c2) -> return $ c1 <= c2
    (Just _,Nothing) -> return True
    (Nothing,Just _) -> return False
    (Nothing,Nothing) ->
      return $ (ml!(k+1)!k)*(ml!(k+1)!k) + (ml!(k+1)!(k+1))*(ml!(k-1)!(k-1)) < (ml!k!k)*(ml!k!k)

-- | ST to compute the Hermite normal form for the matrix in stMatA.
-- The first argument is also subject to the operation that is executed on stMatA, so one gets the transformation unimodular matrix by passing the identity.
-- See the source of hnfLLL function for an example of usage.
hnfLLLST :: (LLLQualified t) => STMatrix s t -> STMatrix s t -> ST s ()
hnfLLLST stMatU stMatA = do
  m <- LA.rows <$> freezeMatrix stMatA
  newLambda <- thawMatrix $ LA.ident (m+1)
  newC1 <- newSTRef Nothing
  newC2 <- newSTRef Nothing
  let resRef = HNFLLL stMatA stMatU newLambda newC1 newC2
  kRef <- newSTRef (1 :: Int)
  whileM_ ((<m) <$> readSTRef kRef) $ do
    k <- readSTRef kRef
    reduceLLL k (k-1) resRef
    shouldSwap <- shouldSwapST k resRef
    if shouldSwap
      then swapLLL k resRef >> when (k>1) (writeSTRef kRef (k-1))
      else (forM_ [0..k-2] $ \i -> reduceLLL k i resRef) >> writeSTRef kRef (k+1)
  setMatrix stMatA 0 0 =<< (LA.flipud <$> freezeMatrix stMatA)
  setMatrix stMatU 0 0 =<< (LA.flipud <$> freezeMatrix stMatU)

-- | Execute the LLL-based algorithm on a given matrix.
hnfLLL :: LLLQualified t => LA.Matrix t -> (LA.Matrix t, LA.Matrix t)
hnfLLL mtx = runST $ do
  stMatA <- thawMatrix mtx
  stMatU <- thawMatrix $ LA.ident (LA.rows mtx)
  hnfLLLST stMatU stMatA
  h <- freezeMatrix stMatA
  u <- freezeMatrix stMatU
  return (u,h)

max_stars_repo_name: Junology/linkedgram

max_issues_repo_name: Junology/linkedgram

max_forks_repo_name: Junology/linkedgram