text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE IncoherentInstances  #-}
{-# LANGUAGE StandaloneDeriving  #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ForeignFunctionInterface #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Data.Vector.Image.IO (
  readImg
, writeImg
, writeBmp
, writePng
, writeJpeg
) where

import Prelude hiding(map)
import qualified Prelude as P
import qualified Data.List as L
import Data.Vector.Storable hiding(forM_,map,replicate,fromList,toList,(!?),(!))
import Data.Maybe
import Numeric.LinearAlgebra.Data hiding (fromList,toList,(!))
import Numeric.LinearAlgebra.HMatrix hiding (fromList,toList,(!))
import qualified Data.Vector.Storable as V
import Foreign
import GHC.Real
import Foreign.C.String
import Foreign.C.Types
import qualified Control.Monad as C
import Data.Vector.Image.Color.RGB
import Data.Vector.Image.Image

type WriteImgFunc = Ptr CChar -- ^ Filename
                 -> Ptr RGB8  -- ^ Image Vector
                 -> CInt      -- ^ Width of Image
                 -> CInt      -- ^ Height of Image
                 -> IO CInt   -- ^ When write fails, this returns negative value. When success, this returns 0.

foreign import ccall unsafe "readImg" c_readImg :: Ptr CChar -> Ptr (Ptr RGB8) -> Ptr CInt -> Ptr CInt -> IO CInt
foreign import ccall unsafe "&freeImg" c_p_freeImg :: FunPtr(Ptr RGB8 -> IO ())
foreign import ccall unsafe "writeBmp" c_writeBmp :: WriteImgFunc
foreign import ccall unsafe "writePng" c_writePng :: WriteImgFunc
foreign import ccall unsafe "writeJpeg" c_writeJpeg :: WriteImgFunc

readImg :: String  -- ^ Filename
        -> IO (Either Int Image8)
readImg file = do
  withCString file $ \ cfile -> do
    alloca $ \cp -> do
      alloca $ \cw -> do
        alloca $ \ch -> do
          r <- fmap fromIntegral $ c_readImg cfile cp cw ch
          w <- fmap fromIntegral $ peek cw
          h <- fmap fromIntegral $ peek ch
          p <- peek cp :: IO (Ptr RGB8)
          let n = w*h
          fp <- newForeignPtr c_p_freeImg p
          return $ if r < 0
                   then Left r
                   else Right $  Image w h (unsafeFromForeignPtr fp 0 n)

writeImg' :: WriteImgFunc -> String -> Image8 -> IO (Either Int ())
writeImg' func file (Image w h dat) = do
  r <- fmap fromIntegral $ withCString file $ \ cfile -> do
    unsafeWith dat $ \p -> do
      func cfile p (fromIntegral w) (fromIntegral h)
  return $ if r < 0 then Left r else Right ()

writeImg :: String   -- ^ Filename
         -> Image8  -- ^ Image Data
         -> IO (Either Int ())
writeImg file img | L.isSuffixOf ".bmp" file = writeBmp file img
                  | L.isSuffixOf ".BMP" file = writeBmp file img
                  | L.isSuffixOf ".png" file = writePng file img
                  | L.isSuffixOf ".Png" file = writePng file img
                  | L.isSuffixOf ".jpg" file = writeJpeg file img
                  | L.isSuffixOf ".JPG" file = writeJpeg file img
                  | otherwise = return $ Left (-1)

writeBmp :: String   -- ^ Filename
         -> Image8  -- ^ Image Data
         -> IO (Either Int ())
writeBmp = writeImg' c_writeBmp

writeJpeg :: String  -- ^ Filename
          -> Image8 -- ^ Image Data
          -> IO (Either Int ())
writeJpeg = writeImg' c_writeJpeg

writePng :: String  -- ^ Filename
         -> Image8 -- ^ Image Data
         ->  IO (Either Int ())
writePng = writeImg' c_writePng

ext: hs

max_stars_repo_name: junjihashimoto/img-vector

max_issues_repo_name: junjihashimoto/img-vector

max_forks_repo_name: junjihashimoto/img-vector


text:
{-|

Copyright:
  This file is part of the package zxcvbn-hs. It is subject to the
  license terms in the LICENSE file found in the top-level directory
  of this distribution and at:

    https://code.devalot.com/sthenauth/zxcvbn-hs

  No part of this package, including this file, may be copied,
  modified, propagated, or distributed except according to the terms
  contained in the LICENSE file.

License: MIT

-}
module Text.Password.Strength.Internal.Math
  ( variations
  , variations'
  , bruteForce
  , caps
  ) where

--------------------------------------------------------------------------------
-- Library Imports:
import Control.Lens ((^.))
import Data.Char (isUpper)
import qualified Data.Text as Text
import Numeric.SpecFunctions (choose)

--------------------------------------------------------------------------------
-- Project Imports:
import Text.Password.Strength.Internal.Token

--------------------------------------------------------------------------------
-- | Equation 2, section 4, page 163 (8/18)
--
-- NOTE: The other implementations don't seem to divide the sum by two
--       but the equation in the paper clearly does.
variations :: Int -> Int -> Integer
variations 0 _ = 1
variations _ 0 = 1
variations u l = max 1 (sum (floor . choose (u+l) <$> [1 .. min u l]) `div` 2)

--------------------------------------------------------------------------------
-- | Like equation 2, but modified for l33t and keyboard variations.
--
-- This equation does not appear in the 2016 paper although it is
-- mentioned.  Therefore it was adapted from the CoffeeScript and
-- Python implementations.
variations' :: Int -> Int -> Integer
variations' 0 _ = 2
variations' _ 0 = 2
variations' u l = variations u l

--------------------------------------------------------------------------------
-- | Calculate the brute force score for text with the given length.
bruteForce :: Int -> Integer
bruteForce = (10 ^)

--------------------------------------------------------------------------------
-- | Score the use of uppercase letters according to the paper.  This
-- is specified in the paragraph preceding equation 2.
caps :: Token -> Integer -> Integer
caps token n =
  let text       = token ^. tokenChars
      upper      = Text.length (Text.filter isUpper text)
      lower      = Text.length text - upper
      allLower   = lower == Text.length text
      allUpper   = lower == 0
      firstUpper = upper == 1 && Text.all isUpper (Text.take 1 text)
      lastUpper  = upper == 1 && Text.all isUpper (Text.takeEnd 1 text)
  in case () of
    () | allLower   -> n
       | firstUpper -> n * 2
       | lastUpper  -> n * 2
       | allUpper   -> n * 2
       | otherwise  -> n * variations upper lower

--------------------------------------------------------------------------------
-- NOTE: Equation 3 is in Keybaord.hs.


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE NamedFieldPuns #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE ConstraintKinds #-}

module Language.Grappa.Interp.ProbFun where

import Data.Vector (Vector)
import qualified Data.Vector as V
import qualified Numeric.Log as Log
import qualified Numeric.LinearAlgebra as M
import Data.Functor.Compose
import Data.Aeson

import Language.Grappa.Distribution
import Language.Grappa.Interp
import Language.Grappa.GrappaInternals
import Language.Grappa.Interp.PVIE
import GHC.Exts (IsList(..))

import Debug.Trace

-- | Type tag for the representation that views @'Dist' a@ as a function from
-- @a@ to 'Prob'
data ProbFunRepr

type family ProbFunReprF a where
  ProbFunReprF (a -> b) = (ProbFunReprF a -> ProbFunReprF b)
  ProbFunReprF (Dist a) = ProbFunReprF a -> Prob
  ProbFunReprF (ADT adt) = adt (GExpr ProbFunRepr) (ADT adt)
  ProbFunReprF (Vector a) = Vector (ProbFunReprF a)
  ProbFunReprF RVector = RVector
  ProbFunReprF RMatrix = RMatrix
  ProbFunReprF ProbVector = ProbVector
  ProbFunReprF ProbMatrix = ProbMatrix
  ProbFunReprF (VIDist a) = (VIDistFamExpr (ProbFunReprF a))
  ProbFunReprF VISize = VIDim
  ProbFunReprF a = a

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

instance ValidRepr ProbFunRepr where
  type GVExprRepr ProbFunRepr a = ProbFunReprF a
  type GStmtRepr ProbFunRepr a = Prob

  interp__'projTupleStmt (GExpr tup) k = k tup

  interp__'vInjTuple tup = GVExpr $ mapADT (GExpr . unGVExpr) tup
  interp__'vProjTuple (GVExpr ve) k = k $ mapADT (GVExpr . unGExpr) ve

  interp__'vwild _ = error "ProbFunRepr: wildcard variables not allowed!"
  interp__'vlift (GExpr e) k = k (GVExpr e)

  interp__'return (GExpr _) = GStmt 1
  interp__'let rhs body = body rhs
  interp__'sample (GExpr d) (GVExpr v) k =
    -- FIXME: evaluate these two sides in parallel!
    GStmt $ d v * unGStmt (k $ GExpr v)

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

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

instance TraversableADT adt => Interp__ADT ProbFunRepr adt where
  interp__'vInjADT adt = GVExpr (mapADT (GExpr . unGVExpr) adt)
  interp__'vProjMatchADT (GVExpr adt) _ matcher k_succ k_fail =
    if applyCtorMatcher matcher adt then
      k_succ (mapADT (GVExpr . unGExpr) adt)
    else k_fail

instance EmbedRepr ProbFunRepr Bool where
  embedRepr = GExpr
instance EmbedRepr ProbFunRepr Int where
  embedRepr = GExpr
instance EmbedRepr ProbFunRepr R where
  embedRepr = GExpr
instance EmbedRepr ProbFunRepr Prob where
  embedRepr = GExpr

instance IsList (ListF a (GExpr ProbFunRepr) (ADT (ListF a))) where
  type Item (ListF a (GExpr ProbFunRepr) (ADT (ListF a))) = ProbFunReprF a
  fromList [] = Nil
  fromList (x:xs) = (Cons (GExpr x) (GExpr (fromList xs)))
  toList Nil = []
  toList (Cons (GExpr x) (GExpr xs)) = x : toList xs

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

instance Interp__'ifThenElse ProbFunRepr where
  interp__'ifThenElse (GExpr c) t e = if c then t else e

instance Interp__'vmatchSwitch ProbFunRepr where
  interp__'vmatchSwitch (GExpr i) stmts = stmts !! i

instance Interp__not ProbFunRepr where
  interp__not = GExpr not

instance Interp__'amp'amp ProbFunRepr where
  interp__'amp'amp = GExpr (&&)

instance Interp__'bar'bar ProbFunRepr where
  interp__'bar'bar = GExpr (||)

instance (Eq a, Eq (ProbFunReprF a)) =>
         Interp__'eq'eq ProbFunRepr a where
  interp__'eq'eq = GExpr (==)

instance (Ord a, Ord (ProbFunReprF a)) =>
         Interp__'lt ProbFunRepr a where
  interp__'lt = GExpr (<)

instance (Ord a, Ord (ProbFunReprF a)) =>
         Interp__'gt ProbFunRepr a where
  interp__'gt = GExpr (>)

instance (Ord a, Ord (ProbFunReprF a)) =>
         Interp__'lt'eq ProbFunRepr a where
  interp__'lt'eq = GExpr (<=)

instance (Ord a, Ord (ProbFunReprF a)) =>
         Interp__'gt'eq ProbFunRepr a where
  interp__'gt'eq = GExpr (>=)

instance (Ord a, Ord (ProbFunReprF a)) =>
         Interp__min ProbFunRepr a where
  interp__min = GExpr min

instance (Ord a, Ord (ProbFunReprF a)) =>
         Interp__max ProbFunRepr a where
  interp__max = GExpr max

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

instance (Num a, Num (ProbFunReprF a)) =>
         Interp__'plus ProbFunRepr a where
  interp__'plus = GExpr (+)

instance (Num a, Num (ProbFunReprF a)) =>
         Interp__'minus ProbFunRepr a where
  interp__'minus = GExpr (-)

instance (Num a, Num (ProbFunReprF a)) =>
         Interp__'times ProbFunRepr a where
  interp__'times = GExpr (*)

instance (Num a, Num (ProbFunReprF a)) =>
         Interp__negate ProbFunRepr a where
  interp__negate = GExpr negate

instance (Num a, Num (ProbFunReprF a)) =>
         Interp__abs ProbFunRepr a where
  interp__abs = GExpr abs

instance (Num a, Num (ProbFunReprF a)) =>
         Interp__signum ProbFunRepr a where
  interp__signum = GExpr signum

instance (Num a, Num (ProbFunReprF a)) =>
         Interp__fromInteger ProbFunRepr a where
  interp__fromInteger = GExpr fromInteger

instance (Num a, Num (ProbFunReprF a)) =>
         Interp__'integer ProbFunRepr a where
  interp__'integer n = GExpr (fromInteger n)

instance (Interp__'integer ProbFunRepr a,
          Eq (ProbFunReprF a))
         => Interp__'eqInteger ProbFunRepr a where
  interp__'eqInteger (GExpr x) (GExpr y) = GExpr (x == y)

instance (Integral a, Integral (ProbFunReprF a)) =>
         Interp__quot ProbFunRepr a where
  interp__quot = GExpr quot

instance (Integral a, Integral (ProbFunReprF a)) =>
         Interp__rem ProbFunRepr a where
  interp__rem = GExpr rem

instance (Integral a, Integral (ProbFunReprF a)) =>
         Interp__div ProbFunRepr a where
  interp__div = GExpr div

instance (Integral a, Integral (ProbFunReprF a)) =>
         Interp__mod ProbFunRepr a where
  interp__mod = GExpr mod

instance (Integral a, Integral (ProbFunReprF a)) =>
         Interp__toInteger ProbFunRepr a where
  interp__toInteger = GExpr toInteger

instance (Fractional a, Fractional (ProbFunReprF a)) =>
         Interp__'div ProbFunRepr a where
  interp__'div = GExpr (/)

instance (Fractional a, Fractional (ProbFunReprF a)) =>
         Interp__recip ProbFunRepr a where
  interp__recip = GExpr recip

instance (Fractional a, Fractional (ProbFunReprF a)) =>
         Interp__fromRational ProbFunRepr a where
  interp__fromRational = GExpr fromRational

instance (Fractional a, Fractional (ProbFunReprF a)) =>
         Interp__'rational ProbFunRepr a where
  interp__'rational n = GExpr (fromRational n)

instance (Interp__'rational ProbFunRepr a,
          Eq (ProbFunReprF a))
         => Interp__'eqRational ProbFunRepr a where
  interp__'eqRational (GExpr x) (GExpr y) = GExpr (x == y)

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__pi ProbFunRepr a where
  interp__pi = GExpr pi

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__exp ProbFunRepr a where
  interp__exp = GExpr exp

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__log ProbFunRepr a where
  interp__log = GExpr log

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__sqrt ProbFunRepr a where
  interp__sqrt = GExpr sqrt

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__'times'times ProbFunRepr a where
  interp__'times'times = GExpr (**)

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__logBase ProbFunRepr a where
  interp__logBase = GExpr logBase

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__sin ProbFunRepr a where
  interp__sin = GExpr sin

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__cos ProbFunRepr a where
  interp__cos = GExpr cos

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__tan ProbFunRepr a where
  interp__tan = GExpr tan

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__asin ProbFunRepr a where
  interp__asin = GExpr asin

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__acos ProbFunRepr a where
  interp__acos = GExpr acos

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__atan ProbFunRepr a where
  interp__atan = GExpr atan

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__sinh ProbFunRepr a where
  interp__sinh = GExpr sinh

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__cosh ProbFunRepr a where
  interp__cosh = GExpr cosh

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__tanh ProbFunRepr a where
  interp__tanh = GExpr tanh

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__asinh ProbFunRepr a where
  interp__asinh = GExpr asinh

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__acosh ProbFunRepr a where
  interp__acosh = GExpr acosh

instance (Floating a, Floating (ProbFunReprF a)) =>
         Interp__atanh ProbFunRepr a where
  interp__atanh = GExpr atanh

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

instance Interp__realToProb ProbFunRepr where
  interp__realToProb = GExpr $ \x ->
    if x < 0 then 0 else rToProb x

instance Interp__logRealToProb ProbFunRepr where
  interp__logRealToProb = GExpr logRToProb

instance Interp__probToReal ProbFunRepr where
  interp__probToReal = GExpr probToR

instance Interp__probToLogReal ProbFunRepr where
  interp__probToLogReal = GExpr probToLogR

instance Interp__gammaProb ProbFunRepr where
  interp__gammaProb = GExpr (Prob . Log.Exp . logGamma)

instance Interp__digamma ProbFunRepr where
  interp__digamma = GExpr digamma

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

instance (Show a, Show (ProbFunReprF a)) =>
         Interp__gtrace ProbFunRepr a b where
  interp__gtrace = GExpr gtrace

instance Interp__gerror ProbFunRepr a where
  interp__gerror = GExpr gerror

----------------------------------------------------------------------
-- Interpreting Distributions
----------------------------------------------------------------------

instance Interp__delta ProbFunRepr Int where
  interp__delta = GExpr $ \i x -> if x == i then 1 else 0

instance Interp__normal ProbFunRepr where
  interp__normal = GExpr normalDensity

instance Interp__uniform ProbFunRepr where
  interp__uniform = GExpr uniformDensity

instance Interp__exponential ProbFunRepr where
  interp__exponential = GExpr $ \rate x -> Prob $ exponentialDensity rate x

instance Interp__gamma ProbFunRepr where
  interp__gamma = GExpr $ \k theta x -> Prob $ gammaDensity k theta x

instance Interp__beta ProbFunRepr where
  interp__beta = GExpr $ \alpha beta x -> Prob $ betaDensity alpha beta x

instance Interp__betaProb ProbFunRepr where
  interp__betaProb = GExpr $ \alpha beta p ->
    Prob $ betaDensityLog alpha beta (fromProb p)

instance Interp__dirichlet ProbFunRepr where
  interp__dirichlet = GExpr $ \alphas xs ->
    Prob $ dirichletDensity (toList alphas) (toList xs)

instance Interp__dirichletProb ProbFunRepr where
  interp__dirichletProb = GExpr $ \alphas xs ->
    Prob $ dirichletDensityLog (toList alphas) (map fromProb $ toList xs)

instance Interp__dirichletV ProbFunRepr where
  interp__dirichletV = GExpr $ \alphas xs ->
    dirichletDensityV alphas xs

instance Interp__dirichletPV ProbFunRepr where
  interp__dirichletPV = GExpr $ \alphas xs ->
    dirichletDensityPV alphas xs

instance Interp__categorical ProbFunRepr where
  interp__categorical = GExpr $ \probs x ->
    if x >= length (toList probs) then
      trace ("Categorical: unexpected value x = " ++ show x) 0
    else
      (toList probs) !! x

instance Interp__poisson ProbFunRepr where
  interp__poisson = GExpr $ \rate k ->
    if rate < 0 then error "Poisson: rate < 0!" else
      if rate == 0 then
        if k == 0 then 1 else 0
      else
        if k < 0 then 0 else
          Prob $ Log.Exp $
          fromIntegral k * log rate - rate - logGamma (fromIntegral $ k + 1)

instance Interp__binary_mixture ProbFunRepr a where
  interp__binary_mixture = GExpr $ \p d1 d2 x ->
    p * d1 x + (1 - p) * d2 x

instance Interp__scaleDist ProbFunRepr a where
  interp__scaleDist = GExpr $ \p d x -> p * d x

instance Interp__ctorDist__ListF ProbFunRepr where
  interp__ctorDist__Nil = GExpr $ \d xs ->
    case xs of
      Nil -> d Tuple0
      _ -> error "Unexpected Cons!"
  interp__ctorDist__Cons = GExpr $ \d xs ->
    case xs of
      Cons x xs' -> d (Tuple2 x xs')
      _ -> error "Unexpected Nil!"

instance Interp__adtDist__ListF ProbFunRepr where
  interp__adtDist__ListF = GExpr $ \probNil dNil probCons dCons xs ->
    case xs of
      Nil -> probNil * dNil Tuple0
      Cons x xs' -> probCons * dCons (Tuple2 x xs')

instance Interp__arbitrary ProbFunRepr a where
  interp__arbitrary = GExpr $ \_ -> 1

----------------------------------------------------------------------
-- Interpreting Vectors
----------------------------------------------------------------------

instance Interp__vec_head ProbFunRepr a where
  interp__vec_head = GExpr V.head

instance Interp__vec_tail ProbFunRepr a where
  interp__vec_tail = GExpr V.tail

instance Interp__vec_nth ProbFunRepr a where
  interp__vec_nth = GExpr (V.!)

instance Interp__vec_length ProbFunRepr a where
  interp__vec_length = GExpr V.length

instance Interp__vec_map ProbFunRepr a b where
  interp__vec_map = GExpr V.map

instance Interp__vec_foldr ProbFunRepr a b where
  interp__vec_foldr = GExpr V.foldr'

instance Interp__vec_iid ProbFunRepr a where
  interp__vec_iid = GExpr $ \n d xs ->
    if V.length xs == n then
      product $ V.map d xs
    else 0

{-
instance Interp__vec_dist ProbFunRepr where
  interp__vec_dist = GExpr $ \d xs -> d (V.toList xs)
-}

instance Interp__vec_nil_dist ProbFunRepr a where
  interp__vec_nil_dist = GExpr $ \xs ->
    if V.null xs then 1 else 0

instance Interp__vec_cons_dist ProbFunRepr a where
  interp__vec_cons_dist = GExpr $ \d xs ->
    if V.null xs then 0 else
      d (Tuple2 (GExpr $ V.head xs) (GExpr $ V.tail xs))

----------------------------------------------------------------------
-- Interpreting unboxed vectors and matrices
----------------------------------------------------------------------

instance Interp__lengthV ProbFunRepr where
  interp__lengthV = GExpr lengthV

instance Interp__atV ProbFunRepr where
  interp__atV = GExpr atV

instance Interp__generateV ProbFunRepr where
  interp__generateV = GExpr generateV

instance Interp__sumV ProbFunRepr where
  interp__sumV = GExpr sumV

instance Interp__rowsM ProbFunRepr where
  interp__rowsM = GExpr rowsM

instance Interp__colsM ProbFunRepr where
  interp__colsM = GExpr colsM

instance Interp__atM ProbFunRepr where
  interp__atM = GExpr atM

instance Interp__fromRowsM ProbFunRepr where
  interp__fromRowsM =
    GExpr (fromRowsM . map unGExpr . toHaskellListF unGExpr)

instance Interp__fromColsM ProbFunRepr where
  interp__fromColsM =
    GExpr (fromColsM . map unGExpr . toHaskellListF unGExpr)

instance Interp__toRowsM ProbFunRepr where
  interp__toRowsM =
    GExpr (fromHaskellListF GExpr . map GExpr . toRowsM)

instance Interp__toColsM ProbFunRepr where
  interp__toColsM =
    GExpr (fromHaskellListF GExpr . map GExpr . toColsM)

instance Interp__mulM ProbFunRepr where
  interp__mulM = GExpr mulM

instance Interp__mulMV ProbFunRepr where
  interp__mulMV = GExpr mulMV

instance Interp__mulVM ProbFunRepr where
  interp__mulVM = GExpr mulVM

----------------------------------------------------------------------
-- Interpreting unboxed vectors and matrices of probabilities
----------------------------------------------------------------------

instance Interp__lengthPV ProbFunRepr where
  interp__lengthPV = GExpr lengthPV

instance Interp__atPV ProbFunRepr where
  interp__atPV = GExpr atPV

instance Interp__generatePV ProbFunRepr where
  interp__generatePV = GExpr generatePV

instance Interp__sumPV ProbFunRepr where
  interp__sumPV = GExpr sumPV

instance Interp__rowsPM ProbFunRepr where
  interp__rowsPM = GExpr rowsPM

instance Interp__colsPM ProbFunRepr where
  interp__colsPM = GExpr colsPM

instance Interp__atPM ProbFunRepr where
  interp__atPM = GExpr atPM

instance Interp__fromRowsPM ProbFunRepr where
  interp__fromRowsPM =
    GExpr (fromRowsPM . map unGExpr . toHaskellListF unGExpr)

instance Interp__fromColsPM ProbFunRepr where
  interp__fromColsPM =
    GExpr (fromColsPM . map unGExpr . toHaskellListF unGExpr)

instance Interp__toRowsPM ProbFunRepr where
  interp__toRowsPM =
    GExpr (fromHaskellListF GExpr . map GExpr . toRowsPM)

instance Interp__toColsPM ProbFunRepr where
  interp__toColsPM =
    GExpr (fromHaskellListF GExpr . map GExpr . toColsPM)

instance Interp__mulPM ProbFunRepr where
  interp__mulPM = GExpr mulPM

instance Interp__mulPMV ProbFunRepr where
  interp__mulPMV = GExpr mulPMV

instance Interp__mulPVM ProbFunRepr where
  interp__mulPVM = GExpr mulPVM

----------------------------------------------------------------------
-- Distributions involving unboxed vectors and matrices of probabilities
----------------------------------------------------------------------

instance Interp__categoricalPV ProbFunRepr where
  interp__categoricalPV = GExpr $ \probs x ->
    if x >= lengthPV probs then
      trace ("Categorical: unexpected value x = " ++ show x) 0
    else
      atPV probs x

instance Interp__categoricalV ProbFunRepr where
  interp__categoricalV = GExpr $ \probs x ->
    if x >= lengthV probs then
      trace ("Categorical: unexpected value x = " ++ show x) 0
    else
      realToProb (atV probs x)

{-
instance Interp__mv_normal ProbFunRepr where
  interp__mvNormal =
-}

----------------------------------------------------------------------
-- Interpreting VI Distribution Families
----------------------------------------------------------------------

instance Interp__withVISize ProbFunRepr where
  interp__withVISize = GExpr $ \f -> bindVIDimFamExpr f

instance Interp__viNormal ProbFunRepr where
  interp__viNormal = GExpr normalVIFamExpr

instance Interp__viUniform ProbFunRepr where
  interp__viUniform = GExpr uniformVIFamExpr

instance Interp__viGamma ProbFunRepr where
  interp__viGamma = GExpr gammaVIFamExpr

instance Interp__viGammaProb ProbFunRepr where
  interp__viGammaProb = GExpr gammaProbVIFamExpr

instance Interp__viBeta ProbFunRepr where
  interp__viBeta = GExpr betaVIFamExpr

instance Interp__viBetaProb ProbFunRepr where
  interp__viBetaProb = GExpr betaProbVIFamExpr

instance Interp__viCategorical ProbFunRepr where
  interp__viCategorical = GExpr categoricalVIFamExpr

instance Interp__viDirichlet ProbFunRepr where
  interp__viDirichlet =
    GExpr $ \sz ->
    xformVIDistFamExpr fromList toList (dirichletVIFamExpr sz)

instance Interp__viDirichletProb ProbFunRepr where
  interp__viDirichletProb =
    GExpr $ \sz ->
    xformVIDistFamExpr fromList toList (dirichletProbVIFamExpr sz)

instance Interp__viDirichletV ProbFunRepr where
  interp__viDirichletV =
    GExpr $ \sz ->
    xformVIDistFamExpr (fromListV . fromList) (toList . toListV)
    (dirichletVIFamExpr sz)

instance Interp__viDirichletPV ProbFunRepr where
  interp__viDirichletPV =
    GExpr $ \sz ->
    xformVIDistFamExpr (fromListPV . fromList) (toList . toListPV)
    (dirichletProbVIFamExpr sz)

instance (Eq (ProbFunReprF a), GrappaShow (ProbFunReprF a)) =>
         Interp__viDelta ProbFunRepr a where
  interp__viDelta = GExpr deltaVIFamExpr

instance (Eq (ProbFunReprF a), GrappaShow (ProbFunReprF a),
          FromJSON (ProbFunReprF a)) =>
         Interp__viJSONInput ProbFunRepr a where
  interp__viJSONInput = GExpr readJSONVIDistFamExpr

instance FromJSON (ProbFunReprF a) =>
         Interp__viMappedJSONInput ProbFunRepr a b where
  interp__viMappedJSONInput = GExpr mapJSONVIDistFamExpr

instance Interp__viTuple0 ProbFunRepr where
  interp__viTuple0 = GExpr $ tupleVIFamExpr Tuple0

instance Interp__viTuple1 ProbFunRepr a where
  interp__viTuple1 = GExpr $ \da ->
    tupleVIFamExpr (Tuple1 $ Compose $ xformVIDistFamExpr GExpr unGExpr da)

instance Interp__viTuple2 ProbFunRepr a b where
  interp__viTuple2 = GExpr $ \da db ->
    tupleVIFamExpr (Tuple2
                    (Compose $ xformVIDistFamExpr GExpr unGExpr da)
                    (Compose $ xformVIDistFamExpr GExpr unGExpr db))

instance Interp__viTuple3 ProbFunRepr a b c where
  interp__viTuple3 = GExpr $ \da db dc ->
    tupleVIFamExpr (Tuple3
                    (Compose $ xformVIDistFamExpr GExpr unGExpr da)
                    (Compose $ xformVIDistFamExpr GExpr unGExpr db)
                    (Compose $ xformVIDistFamExpr GExpr unGExpr dc))

instance Interp__viTuple4 ProbFunRepr a b c d where
  interp__viTuple4 = GExpr $ \da db dc dd ->
    tupleVIFamExpr (Tuple4
                    (Compose $ xformVIDistFamExpr GExpr unGExpr da)
                    (Compose $ xformVIDistFamExpr GExpr unGExpr db)
                    (Compose $ xformVIDistFamExpr GExpr unGExpr dc)
                    (Compose $ xformVIDistFamExpr GExpr unGExpr dd))

instance Interp__viIID ProbFunRepr a where
  interp__viIID =
    GExpr $ \sz d ->
    xformVIDistFamExpr fromList toList (iidVIFamExpr sz d)

instance Interp__viVecDist ProbFunRepr a where
  interp__viVecDist = GExpr vecDistVIFamExpr

instance Interp__viVecIID ProbFunRepr a where
  interp__viVecIID = GExpr vecIIDVIFamExpr

instance Interp__viIIDPV ProbFunRepr where
  interp__viIIDPV = GExpr iidPVVIFamExpr

max_stars_repo_name: kquick/grappa

max_issues_repo_name: kquick/grappa

max_forks_repo_name: kquick/grappa


text:
{-# OPTIONS_GHC  -fno-warn-unused-binds -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-missing-signatures #-}
{-# LANGUAGE FlexibleInstances, ConstraintKinds, ExistentialQuantification, GADTs, RankNTypes, MultiParamTypeClasses   #-}
{-# LANGUAGE ImpredicativeTypes, RankNTypes, UndecidableInstances, FlexibleContexts, TypeSynonymInstances #-}

---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------
-- |
-- | Module : Count-Min Sketch
-- | Creator: Xiao Ling
-- | Created: 12/17/2015
-- |
---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------

module CountMinSketch where

import Control.Monad.Random.Class
import Control.Monad.Random
import Control.Monad.State
import Control.Monad.Reader
import Control.Monad.Identity

import Data.Matrix
import Data.Foldable (toList)
import Data.Hashable
import Data.Conduit
import Data.Sequence (Seq,(|>),update,empty,index)
import qualified Data.Sequence as S
import qualified Data.Conduit.List as Cl

import Core
import Utils
import Statistics

{-----------------------------------------------------------------------------
  Types
------------------------------------------------------------------------------}

-- * The table has width `W` and height `D`
-- * items in input stream comes from universe [n] = {1...n}
-- * and let p = 2^n - 1
-- * hash function i |-> a * i + b mod p mod w
type W        = Int
type D        = Int
type P        = Int
type A        = Int
type B        = Int

-- * Our `Table` is a table of `Counters` of size W x D
type Counters = Seq Int
type Table    = (Counters, W, D)

-- * hash function `HashF` = a * i + b mod p mod w
-- * which maps some `i :: a` to `a`
-- * hash functions `HashFs` is a sequence of hash functions
type HashF    = Int -> Int
type HashFs   = Seq HashF

-- * Our algorithm modifies the `Table`, `HashF`unctions, a `Counter`,
-- * and is equipped with randomness
--type Some m = (MonadState (Table, HashFs) m, MonadState Counter m, MonadRandom m)

type Some m   = (MonadState (Table, HashFs) m, MonadRandom m)

{-----------------------------------------------------------------------------
  Play with the correct data type before jumping into it
------------------------------------------------------------------------------}

-- * declare random a for a_1 ... a_d
as :: MonadRandom m => P -> m [A]
as p = take p <$> getRandomRs (1,p)

{-----------------------------------------------------------------------------
  Subroutines
------------------------------------------------------------------------------}

-- * A vacuous inital condition to be thrown away
em :: (Table,HashFs)
em = (,) (S.replicate 0 0, 0, 0) empty

-- * initalize a table given `eps`ilon and `delt`a
table :: Some m => EpsDelta -> m Table
table (ED eps del) = return $ (S.replicate (w * d) 0, w, d)
    where (w,d) = (ceiling $ 2/eps, ceiling $ log2 (1/del))

-- * pick `a :: A` and `b ::B` of hash function given `p ::P`
-- * that is we pick and b uniformly from [1..p]
ab :: Some m => P -> m (A,B)
ab p = (\(a:b:_) -> (a,b)) <$> take 2 <$> getRandomRs (1,p)

-- * construct a `hashF`unction given `p` and `w`
hashF :: Some m => P -> W -> m HashF
hashF p w = do
    (a,b) <- ab p
    return $ \i -> ((a * i + b) `mod` p) `mod` w

hs :: Some m => m [(A,B)]
hs = replicateM 10 $ ab 100

-- * construct `D` pairwise independent `hashF`unctions
-- * TODO : make sure they're actually different
hashFs :: Some m => D -> P -> W -> m HashFs
hashFs d p = S.replicateM d . hashF p

-- * initalize count min algorithm
inits :: Some m => EpsDelta -> P -> m ()
inits eps p = do
    t@(_,w,d) <- table eps
    hs        <- hashFs d p w
    put (t,hs)

update :: Some m => m ()
update = do
    ((cs,w,d), hs) <- get
    return ()

{-----------------------------------------------------------------------------
  CountMinSketch
------------------------------------------------------------------------------}


text:
{-# LANGUAGE OverloadedLists #-}

module Isumi.Math.FFT.PolynomialSpec
  ( spec
  )
where

import           Isumi.Hspec
import           Test.Hspec

import           Data.Complex
import           Isumi.Math.FFT.Internal.Polynomial

spec :: Spec
spec = do
  describe "modPoly" $ do
    it "returns 5 for (x^3 - 2x^2 - 4) `modPoly` (x - 3)" $ do
      let p = PolyVecRep [1, -2, 0, -4] :: Polynomial Double
      let actual = p `modPoly` PolyExpCoefRep [(1, 1), (0, -3)]
      (actual >>= polyToVecRep) `shouldBeAppx` Just [5]
    it "returns 0 for 1 `modPoly` 1" $ do
      let actual = (PolyVecRep [1] :: Polynomial Double)
            `modPoly` PolyExpCoefRep [(0, 1)]
      (actual >>= polyToVecRep) `shouldBeAppx` Just [1]
    it "returns -2+2i for 4x^3+3x^2+2x+1 `modPoly` x+i" $ do
      let p = PolyVecRep [4, 3, 2, 1] :: Polynomial (Complex Double)
          d = PolyExpCoefRep [(1, 1), (0, 0 :+ 1)]
      (p `modPoly` d) `shouldBeAppx` Just (PolyVecRep [(-2) :+ 2])
  describe "polyMultiply'" $
    it "returns (x^3 + x) for x and (x^2 + 1)" $ do
      let expected = [(3, 1), (1, 1)] :: [(Int, Double)]
      let actual = (1, 1) `polyMultiply'` [(2, 1), (0, 1)]
      expected `shouldBeAppx` actual
  describe "factorDivisor" $ do
    it "factors x^2 - 1 into (x - 1)(x + 1)" $
      factorDivisor (BruunDivisorMinus 2) `shouldBeAppx`
        Just (BruunDivisorMinus 1, BruunDivisorPlus 1 (0 :: Double))
    it "factors x + 1 into Nothing" $
      factorDivisor (BruunDivisorComplex (1 :: Complex Double)) `shouldBeAppx`
        Nothing
    it "factors x^4+1.5x^2+1 into (x^2+(sqrt 0.5)x+1)(x^2-(sqrt 0.5)x+1)" $ do
      let c = sqrt 0.5 :: Double
      factorDivisor (BruunDivisorPlus 4 1.5) `shouldBeAppx`
        Just (BruunDivisorPlus 2 c, BruunDivisorPlus 2 (-c))

max_stars_repo_name: IsumiF/fft-bruun

max_issues_repo_name: IsumiF/fft-bruun

max_forks_repo_name: IsumiF/fft-bruun


text:
import Debug.Trace
-- import Control.Monad.State.Lazy

import System.Environment (getArgs)
import System.IO (readFile)
import Data.Map.Strict (Map, (!), insert, elems, fromList, toList, findWithDefault, size, empty, member, findMin, findMax, singleton, filter)
import qualified Data.Map.Strict as M
--import qualified Data.Array as A
import Data.List (permutations)
import Data.List.Split (splitOn)
import Data.Complex (Complex((:+)), realPart, imagPart)

import UI.NCurses

type Instructions = Map Integer Integer

run :: (Integer, Integer) -> Instructions -> [Integer] -> [Integer]
-- (instructionPointer, base), instructions, input, output
run (i, base) instructions inputs =
  -- traceShow (i, base, instructions) $
  case instr `mod` 100 of
    1 -> run (i+4, base) (insert (addr 3) (arg 1 + arg 2) instructions) inputs
    2 -> run (i+4, base) (insert (addr 3) (arg 1 * arg 2) instructions) inputs
    3 -> run (i+2, base) (insert (addr 1) (head inputs) instructions) $ tail inputs
    4 -> -- traceShow (arg 1) $
      arg 1: run (i+2, base) instructions inputs
    5 -> run (if arg 1 == 0 then i+3 else arg 2, base) instructions inputs
    6 -> run (if arg 1 == 0 then arg 2 else i+3, base) instructions inputs
    7 -> run (i+4, base) (insert (addr 3) (if arg 1 < arg 2 then 1 else 0) instructions) inputs
    8 -> run (i+4, base) (insert (addr 3) (if arg 1 == arg 2 then 1 else 0) instructions) inputs
    9 -> run (i+2, base+arg 1) instructions inputs
    99 -> []
    _ -> error "unknown opcode"
  where instr = instructions!i
        ii x = findWithDefault 0 x instructions
        arg :: Integer -> Integer
        arg n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ ii $ i+n
          1 -> ii $ i+n
          2 -> ii ((ii $ i+n) + base)
          _ -> error "bad argument mode"
        addr n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ i+n
          1 -> error "address in mode 1"
          2 -> --trace "address in mode 2" $
            ii (i+n) + base
          _ -> error $ "address in unknown mode"

instance Ord a => Ord (Complex a) where
  compare a b | ar < br || ar == br && ai < bi = LT
              | ar > br || ar == br && ai > bi = GT
              | otherwise = EQ
    where (ar, ai) = toParts a
          (br, bi) = toParts b

toParts :: Complex a -> (a, a)
toParts c = (realPart c, imagPart c)

toIntParts :: RealFrac a => Complex a -> (Int, Int)
toIntParts c = (round $ realPart c, round $ imagPart c)

fromParts :: (a, a) -> Complex a
fromParts (a,b) = a :+ b

type Game = Map (Integer, Integer) Integer
updateGame :: Game -> [Integer] -> Game
updateGame g [] = g
updateGame g (x:y:t:r) = updateGame (insert (x,y) t g) r
updateGame g _ = error "bad input"

display :: [Integer] -> IO ()
display drawCommands = runCurses $ do
    setEcho False
    w <- defaultWindow
    updateWindow w $ do
        draw drawCommands
    render
    waitFor w (\ev -> ev == EventCharacter 'q' || ev == EventCharacter 'Q')
      where draw [] = do
              moveCursor 0 10
              drawString "Done"
            draw (x:y:t:r) =
              if (x,y) == (-1, 0)
              then do
                moveCursor 0 0
                drawString $ show t
                draw r
              else do
                moveCursor (y+1) x
                drawString (case t of
                              0 -> " "
                              1 -> "W"
                              2 -> "B"
                              3 -> "_"
                              4 -> "*")
                draw r

waitFor :: Window -> (Event -> Bool) -> Curses ()
waitFor w p = loop where
    loop = do
        ev <- getEvent w Nothing
        case ev of
            Nothing -> loop
            Just ev' -> if p ev' then return () else loop

play :: [Integer] -> IO ([Integer])
play drawInstructions = runCurses $ do
  setEcho False
  w <- defaultWindow
  loop w drawInstructions
    where loop :: Window -> [Integer] -> Curses([Integer])
          loop w [] = do
            updateWindow w $ do
              moveCursor 0 10
              drawString "Press q to quit"
            render
            let qloop = do
                  ev <- getEvent w Nothing
                  case ev of
                    Nothing -> qloop
                    Just (EventCharacter 'q') -> return []
                    Just (EventCharacter 'Q') -> return []
                    _ -> qloop
            qloop
          loop w (x:y:t:r) = do
            updateWindow w $ do
              if (x,y) == (-1, 0)
                then do
                moveCursor 0 0
                drawString $ show t
                else do
                moveCursor (y+1) x
                drawString (case t of
                              0 -> " "
                              1 -> "W"
                              2 -> "B"
                              3 -> "_"
                              4 -> "*"
                              _ -> error ("Unknown input " ++ show t)
                           )
            render
            moveEvent <- getEvent w (Just 1)
            let joystick = case moveEvent of
                             Just (EventSpecialKey KeyLeftArrow) -> -1
                             Just (EventSpecialKey KeyRightArrow) -> 1
                             Nothing -> 0
                             Just (EventCharacter 'j') -> -1
                             Just (EventCharacter 'J') -> -1
                             Just (EventCharacter 'k') -> 1
                             Just (EventCharacter 'K') -> 1
                             _ -> 0
            updateWindow w $ do
              moveCursor 0 10
              drawString $ case joystick of
                -1 -> "move left"
                0 -> "no move"
                1 -> "move right"
            render
            -- laterJoysticks <- loop r
            -- return (joystick: laterJoysticks)
            joystick `merge` loop w r
          loop w badQ = error $ "loop: bad q: " ++ show badQ
          merge:: Integer -> Curses([Integer]) -> Curses([Integer])
          merge j ~js' = do
            js <- js'
            return (j:js)

main = do
  [instructionFile] <- getArgs
  instructionStrings <- readFile instructionFile
  let instructions = fromList . zip [0 ..] $ map read $ splitOn "," instructionStrings

  putStrLn "Part 1"
  let moves = run (0,0) instructions []
  putStrLn $ show moves
  let game = updateGame empty moves
  putStrLn $ show ("blocks", size (M.filter (== 2) game))

  -- putStrLn "Part 1b"
  -- display drawCommands
  -- putStrLn "Part 2"
  -- let instructions2 = insert 0 2 instructions -- play for free by setting memory address 0 to 2
  -- let moves = moves >>= (play . run (0,0) instructions2)
  -- -- moves <- play $ run (0,0) instructions2 $ repeat 0
  -- moves' <- moves
  -- putStrLn $ show $ length moves'

  -- let moves = moves >>= (play . run (0,0) instructions2)
  -- moves' <- moves
  -- putStrLn $ show moves'
  -- moves' <- play $ run (0,0) instructions2 []
  -- putStrLn $ show ("moves", moves')
  -- let moves = run (0,0) instructions2 $ repeat 0
  -- putStrLn $ show ("moves", moves)
  -- display moves

  -- let out = run (0,0) instructions $ (startColors!origin) : processOutput
  --     (processOutput, finalp, finald, colors) = process (origin,up, startColors) $ out
  -- putStrLn
  --   $ unlines
  --   $ ["done Part 1"
  --     , show $ ("tiles painted: ", size colors)
  --     , show $ ("finalp", finalp, "finald", robot finald)
  --     ]
  -- putStr $ showColors finalp finald colors

  -- putStrLn "Part 2"
  -- let out = run (0,0) instructions $ (startColors2!origin) : processOutput
  --     (processOutput, finalp, finald, colors) = process (origin,up, startColors2) $ out
  -- putStrLn
  --   $ unlines
  --   $ ["done Part 2"
  --     , show $ ("tiles painted: ", size colors)
  --     , show $ ("finalp", finalp, "finald", robot finald)
  --     ]
  -- putStr $ showColors finalp finald colors


text:
{-# OPTIONS_GHC  -fno-warn-unused-binds -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-missing-signatures #-}
{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, UndecidableInstances, FlexibleContexts, TypeSynonymInstances #-}

---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------
-- |
-- | Module : First attempt at approximate counting
-- | Creator: Xiao Ling
-- | Created: 12/08/2015
-- | TODO   : test standard deviation of alpha, beta, and final version
-- |
---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------

module Morris (
    morris
  , morris'
  , count
  , count'
  ) where

import Data.Random
import Data.Conduit
import Data.List.Split
import qualified Data.Conduit.List as Cl
import Control.Monad.Identity
import System.Environment

import Core
import Statistics

{-----------------------------------------------------------------------------
    I. Morris Algorithm list of counter
------------------------------------------------------------------------------}

-- * Count the number of items in `as` to within `eps` of actual
-- * with confidence `delta`
morris :: Eps -> Delta -> [a] -> IO Counter
morris e d as = runRVar (goMorris e d as) StdRandom

goMorris :: Eps -> Delta -> [a] -> RVar Counter
goMorris e d as = medianOfMeans $ Cl.sourceList as $$ count $ cs
  where cs = replicate (t*m) 0
        t  = round $ 1/(e^2*d)  :: Int
        m  = round . log $ 1/d
        medianOfMeans = fmap median' . (fmap . fmap) mean' . fmap (chunksOf t)

count :: [Counter] -> Sink a RVar [Counter]
count cs = (\c -> 2^(round c) - 1) `ffmap` Cl.foldM (\cs _ -> traverse incr cs) cs
  where ffmap = fmap . fmap

{-----------------------------------------------------------------------------
    II. Morris Algorithm list of list of counter
------------------------------------------------------------------------------}

-- * Count the number of items in `as` to within `eps` of actual
-- * with confidence `delta`
morris' :: Eps -> Delta -> [a] -> IO Counter
morris' e d as = runRVar (goMorris' e d as) StdRandom

-- * 160000
-- * Run on stream inputs `as` for t independent trials for `t = 1/eps^2 * d`,
-- * and `m` times in parallel, for `m = log(1/d)` and take the median
goMorris' :: Eps -> Delta -> [a] -> RVar Counter
goMorris' e d as = medianOfMeans $ Cl.sourceList as $$ count' $ ccs
      where
        medianOfMeans = fmap median' . (fmap . fmap) mean'
        ccs           = replicate m $ replicate t 0
        t             = round $ 1/(e^2*d)
        m             = round . log $ 1/d

-- * Given an m-long list `ccs` of lists (each of which is t-lengthd) of counters,
-- * consume the stream and output result
count' :: [[Counter]] -> Sink a RVar [[Counter]]
count' ccs = (\x -> 2^(round x) - 1) `fffmap` Cl.foldM (\xs _ -> incrs' xs) ccs
  where fffmap = fmap . fmap . fmap

-- * given a list of list of counters toss a coin for each counter and incr
-- * this can be flattened
incrs' :: [[Counter]] -> RVar [[Counter]]
incrs' = sequence . fmap (sequence . fmap incr)

{-----------------------------------------------------------------------------
    III. Utils
------------------------------------------------------------------------------}

-- * Increment a counter `x` with probability 1/2^x
incr :: Counter -> RVar Counter
incr x = do
  h <- toss . coin $ 0.5^(round x)
  return $ if isHead h then (seq () succ x) else seq () x

mean', median' :: (Floating a, Ord a, RealFrac a) => [a] -> Float
mean'   = fromIntegral . round . mean
median' = fromIntegral . round . median

{-----------------------------------------------------------------------------
    IV. Naive Implementation
------------------------------------------------------------------------------}

-- * 2x slower than morris'
-- * Run on stream inputs `xs` for t independent trials for `t = 1/eps`,
-- * and `m` times in parralell where `m = 1/(e^2 * d)`, take the median
morrisNaive :: Eps -> Delta -> [a] -> IO Counter
morrisNaive e d = fmap median' . replicateM m . morrisB t
  where (t,m) = (round $ 1/(e^2*d), round $ 1/d)

-- * Run Morris beta on stream inputs `xs` for `t` independent trials and average
morrisB :: Int -> [a] -> IO Counter
morrisB t =  fmap mean' . replicateM t . morrisA

-- * Run Morris alpha on stream inputs `xs`
morrisA :: [a] -> IO Counter
morrisA xs = flip runRVar StdRandom $ Cl.sourceList xs $$ alpha

-- * A step in morris Algorithm alpha
alpha :: Sink a RVar Counter
alpha = (\x -> 2^(round x) - 1) <$> Cl.foldM (\x _ -> incr x) 0


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

module Trade.Statistics.TradeStatistics where

import qualified Data.Vector as Vec

import Data.Time.Clock (NominalDiffTime)

import qualified Data.List as List

import Data.Ord (comparing)

import qualified Statistics.Sample as Sample

import Trade.Type.DeltaSignal (DeltaSignal(..))
import qualified Trade.Type.DeltaSignal.Algorithm as DSA

import Trade.Type.Yield (LogYield(..), logYield2yield)

import Trade.Report.ToReport (ToReport, toReport)

import qualified Trade.Report.Table as Table

import Trade.Statistics.Statistics (Statistics(..), formatYield, formatStat)

import Trade.Statistics.Algorithm (mean, stdDev)

data TradeStatistics ohlc = TradeStatistics {
  maxPeak :: LogYield ohlc
  , meanPeak :: Statistics NominalDiffTime Double
  , stdDevPeak :: Statistics NominalDiffTime Double
  , maxDrawdown :: LogYield ohlc
  , meanDrawdown :: Statistics NominalDiffTime Double
  , stdDevDrawdown :: Statistics NominalDiffTime Double
  }

tradeStatistics :: [DeltaSignal ohlc] -> TradeStatistics ohlc
tradeStatistics dts =
  let maxs = map DSA.maximum dts
      mins = map DSA.minimum dts
      (maxDts, maxZs) = unzip (map (\(LogYield dt y) -> (dt, y)) maxs)
      maxDtsVec = Vec.map realToFrac (Vec.fromList maxDts)
      maxZsVec = Vec.fromList maxZs
      (minDts, minZs) = unzip (map (\(LogYield dt y) -> (dt, y)) mins)
      minDtsVec = Vec.map realToFrac (Vec.fromList minDts)
      minZsVec = Vec.fromList minZs
  in TradeStatistics {
    maxPeak = List.maximumBy (comparing logYield) maxs
    , meanPeak = Statistics (mean maxDtsVec) (exp (Sample.mean maxZsVec))
    , stdDevPeak = Statistics (stdDev maxDtsVec) (exp (Sample.stdDev maxZsVec))
    , maxDrawdown = List.minimumBy (comparing logYield) mins
    , meanDrawdown =  Statistics (mean minDtsVec) (exp (Sample.mean minZsVec))
    , stdDevDrawdown = Statistics (stdDev minDtsVec) (exp (Sample.stdDev minZsVec))
    }

tradeStatistics2table :: TradeStatistics ohlc -> [[String]]
tradeStatistics2table ts =
  [ ["Trade Statistics", "Yield", "Dur. from trade start"]
  , []
  , "Max. peak" : formatYield (logYield2yield (maxPeak ts))
  , "Mean peak" : formatStat (meanPeak ts)
  , "StdDev peak" : formatStat (stdDevPeak ts)
  , []
  , "Max. drawdown" : formatYield (logYield2yield (maxDrawdown ts))
  , "Mean drawdown" : formatStat (meanDrawdown ts)
  , "StdDev drawdown" : formatStat (stdDevDrawdown ts)
  ]

instance ToReport (TradeStatistics ohlc) where
  toReport = Table.table . tradeStatistics2table

toTradeStatistics ::
  (Functor f) =>
  f [DeltaSignal ohlc] -> f (TradeStatistics ohlc)
toTradeStatistics = fmap tradeStatistics

max_stars_repo_name: fphh/trade

max_issues_repo_name: fphh/trade

max_forks_repo_name: fphh/trade


text:
{-# LANGUAGE DeriveDataTypeable    #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE TypeFamilies          #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Plots.Axis
-- Copyright   :  (C) 2015 Christopher Chalmers
-- License     :  BSD-style (see the file LICENSE)
-- Maintainer  :  Christopher Chalmers
-- Stability   :  experimental
-- Portability :  non-portable
--
-- The 'Axis' is the main data type for "plots". It holds all the
-- necessary infomation to be rendered into a 'Diagram'.
--
----------------------------------------------------------------------------
module Plots.Axis
  ( -- * Axis type
    Axis
  , axes
  , axisPlots
  , currentPlots
  , finalPlots
  , plotModifier
  , axisSize
  , colourBarRange

    -- * Predefined axes
  , r2Axis
  , polarAxis

    -- ** Base space
  , BaseSpace

    -- * Axis plots
  , addPlot
  , addPlotable
  , addPlotable'

    -- * Single axis
  , SingleAxis

    -- * Specific axes
    -- ** x-axis
  , xAxis
  , xLabel
  , xMin
  , xMax

    -- ** y-axis
  , yAxis
  , yLabel
  , yMin
  , yMax

    -- ** r-axis
  , rAxis
  , rLabel
  , rMax

    -- ** theta-axis
  , thetaAxis
  , thetaLabel

    -- ** z-axis
  , zAxis
  , zLabel
  , zMin
  , zMax
  ) where

import           Control.Monad.State
import           Data.Complex
import           Data.Default
import           Data.Typeable

import           Diagrams.Coordinates.Polar
import           Diagrams.Prelude
import           Diagrams.TwoD.Text

import           Plots.Axis.ColourBar
import           Plots.Axis.Grid
import           Plots.Axis.Labels
import           Plots.Axis.Line
import           Plots.Axis.Scale
import           Plots.Axis.Title
import           Plots.Axis.Ticks
import           Plots.Legend
import           Plots.Style
import           Plots.Types

import           Linear

------------------------------------------------------------------------
-- Axis data type
------------------------------------------------------------------------

-- Single axis ---------------------------------------------------------

-- | Render information for a single axis line.
data SingleAxis b v n = SingleAxis
-- note the the v is only present for Style v n
  { saLabel     :: AxisLabel b v n
  , saLine      :: AxisLine v n
  , saTickLabel :: TickLabels b v n
  , saScaling   :: AxisScaling n
  , saGridLines :: GridLines v n
  , saTicks     :: Ticks v n
  , saVisible   :: Bool
  }

type instance V (SingleAxis b v n) = v
type instance N (SingleAxis b v n) = n

instance (TypeableFloat n, Renderable (Text n) b)
    => Default (SingleAxis b V2 n) where
  def = SingleAxis
    { saLabel      = def
    , saLine       = def
    , saTickLabel  = def
    , saGridLines  = def
    , saTicks      = def
    , saScaling    = def
    , saVisible    = True
    }

instance Functor f => HasTicks f (SingleAxis b v n) where
  bothTicks = lens saTicks (\sa ticks -> sa {saTicks = ticks})

instance Functor f => HasMajorTicks f (SingleAxis b v n) where
  majorTicks = bothTicks . majorTicks

instance Functor f => HasMinorTicks f (SingleAxis b v n) where
  minorTicks = bothTicks . minorTicks

instance Functor f => HasAxisLabel f (SingleAxis b v n) b where
  axisLabel = lens saLabel (\sa l -> sa {saLabel = l})

instance Functor f => HasTickLabels f (SingleAxis b v n) b where
  tickLabel = lens saTickLabel (\sa tl -> sa {saTickLabel = tl})

instance Functor f => HasAxisLine f (SingleAxis b v n) where
  axisLine = lens saLine (\sa l -> sa {saLine = l})

instance Functor f => HasGridLines f (SingleAxis b v n) where
  gridLines = lens saGridLines (\sa l -> sa {saGridLines = l})

instance Functor f => HasMajorGridLines f (SingleAxis b v n) where
  majorGridLines = gridLines . majorGridLines

instance Functor f => HasMinorGridLines f (SingleAxis b v n) where
  minorGridLines = gridLines . minorGridLines

instance Functor f => HasAxisScaling f (SingleAxis b v n) where
  axisScaling = lens saScaling (\sa s -> sa {saScaling = s})

instance HasVisibility (SingleAxis b v n) where
  visible = lens saVisible (\sa b -> sa {saVisible = b})

-- singleAxisScale :: Lens' (SingleAxis b v n) AxisScale
-- singleAxisScale = lens saScale (\sa s -> sa {saScale = s})

-- singleAxisBound :: Lens' (SingleAxis b v n) (Bound n)
-- singleAxisBound = lens saBounds (\sa b -> sa {saBounds = b})

------------------------------------------------------------------------
-- Axis type
------------------------------------------------------------------------

-- Base space ----------------------------------------------------------

-- | This family is used so that we can say (Axis Polar) but use V2 for the
--   underlying diagram.
type family BaseSpace (c :: * -> *) :: * -> *

type instance BaseSpace V2      = V2
type instance BaseSpace Complex = V2
type instance BaseSpace Polar   = V2
type instance BaseSpace V3      = V3

-- Axis data type ------------------------------------------------------

-- | Axis is the data type that holds all the necessary information to render
--   a plot. Common 'LensLike's used for the axis (see haddock's
--   instances for a more comprehensive list):
--
--   * 'axisStyle'    - customise the 'AxisStyle'
--   * 'legend'       - customise the 'Legend'
--   * 'colourBar'    - customise the 'ColourBar'
--   * 'currentPlots' - current plots in the 'Axis'
--   * 'finalPlots'   - changes to the plots just before rendering
--   * 'axes'         - changes to each 'SingleAxis'
--
--          * 'xAxis' - the x-axis
--          * 'yAxis' - the y-axis
--          * 'zAxis' - the z-axis
--
--   The following 'LensLike's can be used on the on all the axes by
--   applying it the to 'Axis' or can be used on a 'SingleAxis' by using
--   it in combination with a specific axis (like 'xAxis').
--
--   * 'axisLabel'   - customise the 'MinorTicks'
--   * 'tickLabel'   - customise the 'TickLabels'
--   * 'minorTicks'  - customise the 'MinorTicks'
--   * 'majorTicks'  - customise the 'MajorTicks'
--   * 'gridLines'   - customise the 'GridLines'
--   * 'axisLine'    - customise the 'AxisLine'
--   * 'axisScaling' - customise the 'AxisScaling'
--
--   Plots are usually added to the axis using specific functions for
--   those plots ('Plots.Types.Line.linePlot, 'Plots.Types.Bar.barPlot').
--   These functions use 'addPlotable' to add the plot to the axis.
data Axis b c n = Axis
  { _axisStyle   :: AxisStyle b (BaseSpace c) n
  , _colourBar   :: ColourBar b n
  , _colourBarR  :: (n,n)
  , _legend      :: Legend b n
  , _axisTitle   :: Title b (BaseSpace c) n
  -- , _axisTitle      :: AxisTitle

  , _axisPlots   :: [DynamicPlot b (BaseSpace c) n]
  , _plotModifier :: Endo (StyledPlot b (BaseSpace c) n)

  -- the v in each axis is only used for the style
  , _axes        :: c (SingleAxis b (BaseSpace c) n)
  } deriving Typeable

-- | Lens onto the separate axes of an axis. Allows changing the
--   coordinate system as long as the 'BaseSpace' is the same.
--
-- @
-- 'axes' :: 'Lens'' ('Axis' b c n) (c ('SingleAxis' b v n))
-- @
axes :: (v ~ BaseSpace c, v ~ BaseSpace c')
     => Lens (Axis b c  n)
             (Axis b c' n)
             (c  (SingleAxis b v n))
             (c' (SingleAxis b v n))
axes = lens _axes (\(Axis a1 a2 a3 a4 a5 a6 a7 _) a8 -> Axis a1 a2 a3 a4 a5 a6 a7 a8)

-- | The list of plots currently in the axis.
axisPlots :: BaseSpace c ~ v => Lens' (Axis b c n) [DynamicPlot b v n]
axisPlots = lens _axisPlots (\a ps -> a {_axisPlots = ps})

-- | Traversal over the current plots in the axis.
--
--   For example, to make all 'ScatterPlot's currently in the axis use a
--   'connectingLine', you can write
--
-- @
-- 'finalPlots' . 'connectingLine' .= 'True'
-- @
currentPlots :: BaseSpace c ~ v => Traversal' (Axis b c n) (DynamicPlot b v n)
currentPlots = axisPlots . traversed

-- | Setter over the final plot before the axis is rendered.
--
--   For example, to make all 'ScatterPlot's in the axis use a
--   'connectingLine' (both currently in the axis and ones added later),
--   you can add
--
-- @
-- 'finalPlots' . 'connectingLine' .= 'True'
-- @
--
finalPlots :: BaseSpace c ~ v => Setter' (Axis b c n) (StyledPlot b v n)
finalPlots = sets $ \f a -> a {_plotModifier = _plotModifier a <> Endo f}

-- | Lens onto the modifier set by 'finalPlots'. This gets applied to
--   all plots in the axis, just before they are rendered.
plotModifier :: BaseSpace c ~ v => Lens' (Axis b c n) (Endo (StyledPlot b v n))
plotModifier = lens _plotModifier (\a f -> a {_plotModifier = f})

-- Axis instances ------------------------------------------------------

type instance V (Axis b v n) = BaseSpace v
type instance N (Axis b v n) = n

instance (Applicative f, Traversable c) => HasTicks f (Axis b c n) where
  bothTicks = axes . traverse . bothTicks

instance (Applicative f, Traversable c) => HasMajorTicks f (Axis b c n) where
  majorTicks = axes . traverse . majorTicks

instance (Applicative f, Traversable c) => HasMinorTicks f (Axis b c n) where
  minorTicks = axes . traverse . minorTicks

instance (Applicative f, Traversable c) => HasGridLines f (Axis b c n) where
  gridLines = axes . traverse . gridLines

instance (Applicative f, Traversable c) => HasMajorGridLines f (Axis b c n) where
  majorGridLines = axes . traverse . majorGridLines

instance (Applicative f, Traversable c) => HasMinorGridLines f (Axis b c n) where
  minorGridLines = axes . traverse . minorGridLines

instance (Applicative f, Traversable c) => HasAxisLabel f (Axis b c n) b where
  axisLabel = axes . traverse . axisLabel

instance (Applicative f, Traversable c) => HasTickLabels f (Axis b c n) b where
  tickLabel = axes . traverse . tickLabel

instance (Applicative f, Traversable c) => HasAxisScaling f (Axis b c n) where
  axisScaling = axes . traverse . axisScaling

instance Settable f => HasPlotOptions f (Axis b c n) b where
  plotOptions = finalPlots . plotOptions

instance Settable f => HasPlotStyle f (Axis b c n) b where
  plotStyle = finalPlots . plotStyle

instance HasLegend (Axis b c n) b where
  legend = lens _legend (\a l -> a {_legend = l})

instance HasTitle (Axis b c n) b where
  title = lens _axisTitle (\a t -> a {_axisTitle = t})

-- | The size used for the rendered axis.
axisSize :: (HasLinearMap c, Num n, Ord n) => Lens' (Axis b c n) (SizeSpec c n)
axisSize = axes . column renderSize . iso mkSizeSpec getSpec -- column axisScaling . asSizeSpec -- iso mkSizeSpec getSpec

-- | The range used for the colour bar limits. This is automatically set
--   when using 'heatMap' or 'heatMap''
colourBarRange :: Lens' (Axis b v n) (n,n)
colourBarRange = lens _colourBarR (\a r -> a {_colourBarR = r})

instance HasAxisStyle (Axis b v n) b where
  axisStyle = lens _axisStyle (\a sty -> a {_axisStyle = sty})

instance HasColourBar (Axis b v n) b where
  colourBar = lens _colourBar (\a cb -> a {_colourBar = cb})

-- Axis functions ------------------------------------------------------

-- $plotable
-- The 'Plotable' class defines ways of converting the data type to a
-- diagram for some axis. There are several variants for adding an axis
-- with constraints @('InSpace' v n a, 'Plotable' a b)@:
--
-- @
-- 'addPlotable'  ::           a -> 'PlotState' a b -> 'AxisState' b v n
-- 'addPlotable''   ::           a ->                  'AxisState' b v n
-- @
--
-- The last argument is a 'PlotState' so you can use @do@ notation to
-- make adjustments to the plot. The @L@ suffix stands for \"legend\",
-- it is equivalent of using 'addLegendEntry' in the 'PlotState'. Since
-- legend entries are so common it has it's own suffix. The following
-- are equivalent:
--
-- @
-- myaxis = 'r2Axis' &~ do
--   'addPlotable'' myplot $ do
--     'addLegendEntry' "my plot"
-- @
--
-- Most of the time you won't use these functions directly. However,
-- other plotting functions follow this naming convention where instead
-- of @a@, it takes the data needed to make the plot.

-- | Add a 'Plotable' 'Plot' to an 'Axis'.
addPlot
  :: (InSpace (BaseSpace c) n p, MonadState (Axis b c n) m, Plotable p b)
  => Plot p b -- ^ the plot
  -> m ()     -- ^ add plot to the 'Axis'
addPlot p = axisPlots <>= [DynamicPlot p]

-- | Add something 'Plotable' to the 'Axis' with a stateful modification
--   of the 'Plot'.
addPlotable
  :: (InSpace (BaseSpace c) n p, MonadState (Axis b c n) m, Plotable p b)
  => p -- ^ the raw plot
  -> State (Plot p b) () -- ^ changes to the plot
  -> m () -- ^ add plot to the 'Axis'
addPlotable p s = addPlot $ execState s (mkPlot p)

-- | Simple version of 'AddPlotable' without any changes 'Plot'.
addPlotable'
  :: (InSpace (BaseSpace v) n p, MonadState (Axis b v n) m, Plotable p b)
  => p    -- ^ the raw plot
  -> m () -- ^ add plot to the 'Axis'
addPlotable' p = addPlotable p (return ())

------------------------------------------------------------------------
-- Predefined axes
------------------------------------------------------------------------

-- | The default axis for plots in the 'V2' coordinate system.
r2Axis
  :: (TypeableFloat n,
     Renderable (Text n) b,
     Renderable (Path V2 n) b)
  => Axis b V2 n
r2Axis = Axis
  { _axisStyle  = fadedColours
  , _colourBar  = defColourBar
  , _colourBarR = (0,1)
  , _axisTitle  = def

  , _legend       = def
  , _axisPlots    = []
  , _plotModifier = mempty

  , _axes = pure def
  }

-- The x-axis ----------------------------------------------------------

-- | Lens onto the x-axis of an 'Axis'.
xAxis :: R1 c => Lens' (Axis b c n) (SingleAxis b (BaseSpace c) n)
xAxis = axes . _x

-- | The label for the x-axis. Shorthand for @'xAxis' . 'axisLabelText'@.
xLabel :: R1 c => Lens' (Axis b c n) String
xLabel = xAxis . axisLabelText

-- | The minimum x value for the axis. If the value if 'Nothing' (the
--   'Default'), the bounds will be infered by the plots in the axis.
xMin :: R1 c => Lens' (Axis b c n) (Maybe n)
xMin = xAxis . boundMin

-- | The minimum x value for the axis. If the value if 'Nothing' (the
--   'Default'), the bounds will be infered by the plots in the axis.
xMax :: R1 c => Lens' (Axis b c n) (Maybe n)
xMax = xAxis . boundMax

-- The y-axis ----------------------------------------------------------

-- | Lens onto the y-axis of an 'Axis'.
yAxis :: R2 c => Lens' (Axis b c n) (SingleAxis b (BaseSpace c) n)
yAxis = axes . _y

-- | The label for the y-axis. Shorthand for @'yAxis' . 'axisLabelText'@.
yLabel :: R2 c => Lens' (Axis b c n) String
yLabel = yAxis . axisLabelText

-- | The minimum y value for the axis. If the value if 'Nothing' (the
--   'Default'), the bounds will be infered by the plots in the axis.
yMin :: R2 c => Lens' (Axis b c n) (Maybe n)
yMin = yAxis . boundMin

-- | The minimum y value for the axis. If the value if 'Nothing' (the
--   'Default'), the bounds will be infered by the plots in the axis.
yMax :: R2 c => Lens' (Axis b c n) (Maybe n)
yMax = yAxis . boundMax

-- The z-axis ----------------------------------------------------------

-- | Lens onto the z-axis of an 'Axis'.
zAxis :: R3 c => Lens' (Axis b c n) (SingleAxis b (BaseSpace c) n)
zAxis = axes . _z

-- | The label for the z-axis. Shorthand for @'zAxis' . 'axisLabelText'@.
zLabel :: R3 c => Lens' (Axis b c n) String
zLabel = zAxis . axisLabelText

-- | The minimum z value for the axis. If the value if 'Nothing' (the
--   'Default'), the bounds will be infered by the plots in the axis.
zMin :: R3 c => Lens' (Axis b c n) (Maybe n)
zMin = zAxis . boundMin

-- | The minimum z value for the axis. If the value if 'Nothing' (the
--   'Default'), the bounds will be infered by the plots in the axis.
zMax :: R3 c => Lens' (Axis b c n) (Maybe n)
zMax = zAxis . boundMax

-- The r-axis ----------------------------------------------------------

-- | Lens onto the radial axis of an 'Axis'.
rAxis :: Radial c => Lens' (Axis b c n) (SingleAxis b (BaseSpace c) n)
rAxis = axes . _radial

-- | The label for the radial axis. Shorthand for @'rAxis' . 'axisLabelText'@.
rLabel :: Radial c => Lens' (Axis b c n) String
rLabel = rAxis . axisLabelText

-- | The minimum z value for the axis. If the value if 'Nothing' (the
--   'Default'), the bounds will be infered by the plots in the axis.
-- rMin :: R3 c => Lens' (Axis b c n) (Maybe n)
-- rMin = zAxis . boundMin

-- | The minimum radial value for the axis. If the value if 'Nothing'
--   (the 'Default'), the bounds will be infered by the plots in the
--   axis.
rMax :: Radial c => Lens' (Axis b c n) (Maybe n)
rMax = rAxis . boundMax

-- The theta-axis ------------------------------------------------------

-- | Lens onto the radial axis of an 'Axis'.
thetaAxis :: Circle c => Lens' (Axis b c n) (SingleAxis b (BaseSpace c) n)
thetaAxis = axes . el etheta

-- | The label for the radial axis. Shorthand for @'rAxis' . 'axisLabelText'@.
thetaLabel :: Circle c => Lens' (Axis b c n) String
thetaLabel = thetaAxis . axisLabelText

-- R3 Axis

-- instance (TypeableFloat n, Enum n, Renderable (Text n) b, Renderable (Path V2 n) b)
--     => Default (Axis b V3 n) where
--   def = Axis
--           { _axisTitle      = Nothing
--           , _axisSize       = mkWidth 300
--           , _axisPlots      = []
--           , _axisLegend     = def
--           , _axisTheme      = coolTheme
--           , _axisLinearMap  = isometricProjection
--           , _axisAxisBounds = Bounds $ pure def
--           , _axisGridLines  = pure def
--           , _axisLabels     = pure def
--           , _axisScaling    = pure def
--           , _axisTickLabels = pure def
--           , _axisTicks      = pure def
--           , _axisLines      = pure def
--           }

-- R3 Axis

polarAxis
  :: (TypeableFloat n,
      Renderable (Text n) b,
      Renderable (Path V2 n) b)
  => Axis b Polar n
polarAxis = Axis
  { _axisStyle  = fadedColours
  , _colourBar  = defColourBar
  , _colourBarR = (0,1)
  , _axisTitle  = def

  , _legend       = def
  , _axisPlots    = []
  , _plotModifier = mempty

  , _axes = pure def
  }

max_stars_repo_name: adamConnerSax/plots

max_issues_repo_name: adamConnerSax/plots

max_forks_repo_name: adamConnerSax/plots


text:
{-|
Module      : Language.Rust.Inline.Context
Description : Defines contexts (rules mapping Rust types to Haskell types)
Copyright   : (c) Alec Theriault, 2017
License     : BSD-style
Maintainer  : alec.theriault@gmail.com
Stability   : experimental
Portability : GHC
-}
{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MagicHash #-}

module Language.Fortran.Inline.Context where

import Language.Fortran.Inline.Pretty ( renderType )

import Language.Fortran.Quote         ( ty )
import Language.Fortran.Syntax        ( Ty(..), Abi(..), FnDecl(..),
                                     Arg(..), Mutability(..), Unsafety(..) )

import Language.Haskell.TH

--import Data.Semigroup             hiding (Arg)
import Data.Monoid
import Data.Typeable               ( Typeable )
import Control.Monad               ( void, liftM2 )
import Data.Traversable            ( for )
import Data.List                   ( intercalate )

--import Data.Int                    ( Int8, Int16, Int32, Int64 )
--import Data.Word                   ( Word8, Word16, Word32, Word64 )
import Foreign.C.Types             -- pretty much every type here is used

import GHC.Exts                    ( Char#, Int#, Word#, Float#, Double#,
--                                   ByteArray#
                                   )

import qualified Control.Monad.Fail as Fail
import Eigen.Internal -- CComplex
import qualified Data.Vector.Storable.Mutable as VM
import qualified Data.Vector.Storable as V
--import qualified Debug.Trace as D

import Language.Rust.Data.Position
import Foreign hiding (void)

--import Data.Complex

instance Fail.MonadFail First where
  fail x = error $ "MonadFail First error: " ++ x

--start for Vector Context
type CArray = Ptr

-- Easier on the eyes
type RType = Ty ()
type HType = Type

-- | Represents a prioritized set of rules for mapping Haskell types into Rust
-- ones and vice versa.
--
-- The 'Context' argument encodes the fact that we may need look
-- recursively into the 'Context' again before possibly producing a Haskell
-- type.
newtype Context =
    Context ( [ RType -> Context -> First (Q HType, Maybe (Q RType),String) ]
            -- Given a Rust type in a quasiquote, we need to look up the
            -- corresponding Haskell type (for the FFI import) as well as the
            -- C-compatible Rust type (if the initial Rust type isn't already
            -- @#[repr(C)]@.

            , [ HType -> Context -> First (Q RType) ]
            -- Given a field in a Haskell ADT, we need to figure out which
            -- (not-necessarily @#[repr(C)]@) Rust type normally maps into this
            -- Haskell type.

            , [ String ]
            -- Source for the trait impls of @MarshalTo@
            , String
            )
  deriving (Semigroup, Monoid, Typeable)

-- | Applicative lifting of the 'Context' instance
instance Semigroup (Q Context) where
  (<>) = liftM2 (<>)

-- | Applicative lifting of the 'Context' instance
instance Monoid (Q Context) where
  mappend = (<>)
  mempty = pure mempty

-- | Search in a 'Context' for the Haskell type corresponding to a Rust type.
-- If the Rust type is not C-compatible, also return a C compatible type. It is
-- expected that:
--
--   1. The Haskell type have a 'Storable' instance
--   2. The C-compatible Rust type have the same layout
--
lookupRTypeInContext :: RType -> Context -> First (Q HType, Maybe (Q RType),String)
lookupRTypeInContext rustType context@(Context (rules, _, _,_)) =
  foldMap (\fits -> fits rustType context) rules

isFProcedurePtr :: RType -> Bool
isFProcedurePtr (FProcedurePtr _ _ _ _) = True
isFProcedurePtr _ = False

-- | Search in a 'Context' for the Rust type corresponding to a Haskell type.
-- Looking up the Rust type using 'lookupRTypeInContext' should yield the
-- initial Haskell type again.
lookupHTypeInContext :: HType -> Context -> First (Q RType)
lookupHTypeInContext haskType context@(Context (_, rules, _,_)) =
  foldMap (\fits -> fits haskType context) rules

-- | Partial version of 'lookupRTypeInContext' that fails with an error message
-- if the type is not convertible.
getRTypeInContext :: RType -> Context -> (Q HType, Maybe (Q RType),String)
getRTypeInContext rustType context =
  let x = takeIdContext context
   in case getFirst (lookupRTypeInContext rustType context) of
    Just found -> found
    Nothing -> ( fail $ unwords [ "Could not find information about"
                                , renderType rustType
                                , "in the context " ++ x ++ " : " ++ show rustType
                                ]
               , Nothing, "failed" )

takeIdContext :: Context -> String
takeIdContext (Context (_,_,_,i)) = i

-- | Partial version of 'lookupHTypeInContext' that fails with an error message
-- if the type is not convertible.
getHTypeInContext :: HType -> Context -> Q RType
getHTypeInContext haskType context =
  case getFirst (lookupHTypeInContext haskType context) of
    Just found -> found
    Nothing -> fail $ unwords [ "Could not find information about"
                              , pprint haskType
                              , "in the context"
                              ]

-- | Make a 'Context' consisting of rules to map the Rust types on the left to
-- the Haskell types on the right. The Rust types should all be @#[repr(C)]@
-- and the Haskell types should all be 'Storable'.
mkContext :: String -> [(Ty a, Q HType, Bool,String)] -> Q Context
mkContext idContext tys = do
    tys' <- traverse (\(rt,qht,mkImpl,idSubcon) -> do { ht <- qht; pure (void rt,ht,mkImpl,idSubcon) }) tys
    pure (Context ( map fits tys'
                  , map rev tys'
                  , map impl tys'
                  , idContext
                  ))
  where
    fits (rts, hts, _, i) rt _ | rt == rts = pure (pure hts, Nothing,i)
                            | otherwise = mempty

    rev (rts, hts, _,_) ht _  | ht == hts = pure (pure rts)
                            | otherwise = mempty

    impl (rts, _, mkImpl,_)   | mkImpl = implMarshalInto rts
                            | otherwise = mempty

-- | Make a default @MarshalInto@ trait impl. (An identity impl)
implMarshalInto :: Ty () -> String
implMarshalInto t = unlines [ "impl MarshalInto<" ++ tyStr ++ "> for " ++ tyStr ++ " {"
                             , "  fn marshal(self) -> " ++ tyStr ++ " { self }"
                             , "}"
                             ]
  where tyStr = renderType t

-- | Make a singleton 'Context' consisting of a rule to map the given Rust type
-- to the given Haskell type.
singleton :: Ty a -> Q HType -> Q Context
singleton rts qht = mkContext "singleton" [(rts, qht, True,"singleton")]

-- * Some handy contexts

-- | Types defined in 'Foreign.C.Types' and the 'libc' crate.
--
-- There should be no conversion required here - these have /identical/ memory
-- layouts (since they both promise to have the same memory layout as C) and are
-- passed on the stack.
libc :: Q Context
libc = mkContext "libc" []
--h--  [ ([ty| libc::c_char      |], [t| CChar      |], False) -- char
--h--  , ([ty| libc::c_schar     |], [t| CSChar     |], False) -- signed char
--h--  , ([ty| libc::c_uchar     |], [t| CUChar     |], False) -- unsigned char
--h--  , ([ty| libc::c_short     |], [t| CShort     |], False) -- short
--h--  , ([ty| libc::c_ushort    |], [t| CUShort    |], False) -- unsigned short
--h--  , ([ty| libc::c_int       |], [t| CInt       |], False) -- int
--h--  , ([ty| libc::c_uint      |], [t| CUInt      |], False) -- unsigned int
--h--  , ([ty| libc::c_long      |], [t| CLong      |], False) -- long
--h--  , ([ty| libc::c_ulong     |], [t| CULong     |], False) -- unsigned long
--h--  , ([ty| libc::ptrdiff_t   |], [t| CPtrdiff   |], False) -- ptrdiff_t
--h--  , ([ty| libc::size_t      |], [t| CSize      |], False) -- size_t
--h--  , ([ty| libc::wchar_t     |], [t| CWchar     |], False) -- wchar_t
--h--  , ([ty| libc::c_longlong  |], [t| CLLong     |], False) -- long long
--h--  , ([ty| libc::c_ulonglong |], [t| CULLong    |], False) -- unsigned long long
--h--  , ([ty| libc::boolean_t   |], [t| CBool      |], False) -- bool
--h--  , ([ty| libc::intptr_t    |], [t| CIntPtr    |], False) -- intptr_t
--h--  , ([ty| libc::uintptr_t   |], [t| CUIntPtr   |], False) -- uintptr_t
--h--  , ([ty| libc::intmax_t    |], [t| CIntMax    |], False) -- intmax_t
--h--  , ([ty| libc::uintmax_t   |], [t| CUIntMax   |], False) -- unsigned intmax_t
--h--  , ([ty| libc::clock_t     |], [t| CClock     |], False) -- clock_t
--h--  , ([ty| libc::time_t      |], [t| CTime      |], False) -- time_t
--h--  , ([ty| libc::useconds_t  |], [t| CUSeconds  |], False) -- useconds_t
--h--  , ([ty| libc::suseconds_t |], [t| CSUSeconds |], False) -- suseconds_t
--h--  , ([ty| libc::c_float     |], [t| CFloat     |], False) -- float
--h--  , ([ty| libc::c_double    |], [t| CDouble    |], False) -- double
--h--  , ([ty| libc::FILE        |], [t| CFile      |], True)  -- FILE
--h--  , ([ty| libc::fpos_t      |], [t| CFpos      |], True)  -- fpos_t
--h--  , ([ty| libc::int8_t      |], [t| Int8       |], False) -- int8_t
--h--  , ([ty| libc::int16_t     |], [t| Int16      |], False) -- int16_t
--h--  , ([ty| libc::int32_t     |], [t| Int32      |], False) -- int32_t
--h--  , ([ty| libc::int64_t     |], [t| Int64      |], False) -- int64_t
--h--  , ([ty| libc::uint8_t     |], [t| Word8      |], False) -- uint8_t
--h--  , ([ty| libc::uint16_t    |], [t| Word16     |], False) -- uint16_t
--h--  , ([ty| libc::uint32_t    |], [t| Word32     |], False) -- uint32_t
--h--  , ([ty| libc::uint64_t    |], [t| Word64     |], False) -- uint64_t
--h--  ]

vecCtx :: Q Context
vecCtx = mkContext "vecCtx"
  [ ([ty| vecptr |], [t| VM.IOVector |], True ,"vecptr")
  , ([ty| veclen |], [t| Int |], True         ,"veclen")
  ]
-- | Basic numeric (and similar) Haskell and Rust types.
--
-- There should be no conversion required here as these should have identical
-- memory layouts.
basic :: Q Context
basic = do
  runIO $ putStrLn $  "=====debug:" ++ show [ty| character |]
  mkContext "basic"
    [ ([ty| char  |], [t| Char    |], True, "char ") -- 4 bytes
    , ([ty| i8    |], [t| Int8    |], True, "i8   ")
    , ([ty| i16   |], [t| Int16   |], True, "i16  ")
    , ([ty| i32   |], [t| Int32   |], True, "i32  ")
    , ([ty| i64   |], [t| Int64   |], True, "i64  ")
    , ([ty| u8    |], [t| Word8   |], True, "u8   ")
    , ([ty| u16   |], [t| Word16  |], True, "u16  ")
    , ([ty| u32   |], [t| Word32  |], True, "u32  ")
    , ([ty| u64   |], [t| Word64  |], True, "u64  ")
    , ([ty| f32   |], [t| Float   |], True, "f32  ")
    , ([ty| f64   |], [t| Double  |], True, "f64  ")
    , ([ty| isize |], [t| Int     |], True, "isize")
    , ([ty| usize |], [t| Word    |], True, "usize")
    , ([ty| bool  |], [t| Word8   |], True, "bool ")
    , ([ty| ()    |], [t| ()      |], True, "()   ")
--   Fortran
    , ([ty| integer       |], [t| Int            |], True , "integer     ")
    , ([ty| logical       |], [t| Bool            |], True , "logical     ")
    , ([ty| real          |], [t| Float           |], True , "real        ")
--  , ([ty| complex       |], [t| CComplex Float  |], True , "complex     ")
    , ([ty| complex(kind=8) |], [t| Ptr (CComplex Float) |], True , "complexFloat")
    , ([ty| complex(kind=16) |], [t| Ptr (CComplex Double) |], True , "complexDouble ")
    , ([ty| character     |], [t| Char           |], True , "character   ")
    , ([ty| real(kind=8)  |], [t| Double         |], True , "real(kind=8)")
    , ((FString (Span NoPosition NoPosition)), [t|CChar |], True, "FString")
    {-
    , ([ty| logical(kind=1) |], [t| Int8            |])
    , ([ty| character(len=1)|], [t| CChar           |])
--    , ([ty| integer         |], [t| CInt            |])
--    , ([ty| integer(kind=2) |], [t| CShort          |])
--    , ([ty| integer(kind=4) |], [t| CInt            |])
--    , ([ty| integer(kind=4) |], [t| CLong           |])
--    , ([ty| integer(kind=8) |], [t| CLLong          |])
--    , ([ty| integer(kind=1) |], [t| CChar           |])
--    , ([ty| integer(kind=4) |], [t| CSize           |])
    , ([ty| integer(kind=1) |], [t| Int8            |])
    , ([ty| integer(kind=2) |], [t| Int16           |])
    , ([ty| integer(kind=4) |], [t| Int32           |])
    , ([ty| integer(kind=8) |], [t| Int64           |])
    , ([ty| real(kind=4)    |], [t| Float           |])
    , ([ty| complex(kind=4) |], [t| CComplex Float  |])
    , ([ty| complex(kind=8) |], [t| CComplex Double |])
    -}
    ]

-- | Basic unboxed Haskell types
--
-- TODO: MutableByteArray#
ghcUnboxed :: Q Context
ghcUnboxed = do
  mkContext "ghcUnboxed"
    [ ([ty| char      |], [t| Char#      |], False, "char  ")
    , ([ty| isize     |], [t| Int#       |], False, "isize ")
    , ([ty| usize     |], [t| Word#      |], False, "usize ")
    , ([ty| f32       |], [t| Float#     |], False, "f32   ")
    , ([ty| f64       |], [t| Double#    |], False, "f64   ")
--    , ([ty| *const i8 |], [t| ByteArray# |], False)
    ]

-- | Haskell pointers map onto Rust pointers. Note that unlike Rust, Haskell
-- doesn't really distinguish between pointers pointing to immutable memory from
-- those pointing to to mutable memory, so it is up to the user to enforce this.
--
-- NOTE: pointers will not support pointed types that require an intermediate
--       Rust type.
  {-
fImmutableVectors :: Q Context
fImmutableVectors = do
  vecConT <- [t| V.Vector |]
  pure (Context ([rule], [rev vecConT], [], "fImmutableVectors"))
  where
  rule vec context = do
    FArray _ t _  <- pure vec
    (t', Nothing) <- lookupRTypeInContext t context
    pure ([t| $t' |], Nothing)

  rev vecConT pt context = do
    AppT vecCon t <- pure pt
    if vecCon /= vecConT
      then mempty
      else do
        t' <- lookupHTypeInContext t context
        pure (FArray (-1) <$> t' <*> pure ())
        --pure (Ptr Mutable <$> t' <*> pure ())

fComplex :: Q Context
fComplex = do
  vecConT <- [t| Complex |]
  pure (Context ([rule], [rev vecConT], [], "fComplex"))
  where
  rule vec context = do
    ff <- pure vec
    case ff of
      (FArray _ t _) -> do
        (t', Nothing,i) <- lookupRTypeInContext t context
        pure ([t| Ptr $t' |], Nothing,i)
      _ -> mempty

  rev vecConT pt context = do
    AppT vecCon t <- pure pt
    if vecCon /= vecConT
      then mempty
      else do
        t' <- lookupHTypeInContext t context
        pure (FArray (-1) <$> t' <*> pure ())
        --pure (Ptr Mutable <$> t' <*> pure ())
-}

fVectors :: Q Context
fVectors = do
  vecConT <- [t| VM.MVector |]
  pure (Context ([rule], [rev vecConT], [], "fVectors"))
  where
  rule vec context = do
    ff <- pure vec
    case ff of
      (FArray _ t _) -> do
        (t', Nothing,i) <- lookupRTypeInContext t context
        pure ([t| Ptr $t' |], Nothing,i)
      _ -> mempty

  rev vecConT pt context = do
    AppT vecCon t <- pure pt
    if vecCon /= vecConT
      then mempty
      else do
        t' <- lookupHTypeInContext t context
        pure (FArray (-1) <$> t' <*> pure ())
        --pure (Ptr Mutable <$> t' <*> pure ())

vectors :: Q Context
vectors = do
  vecConT <- [t| V.MVector |]
  pure (Context ([rule], [rev vecConT], [], "vectors"))
  where
  rule vec context = do
    Array t _ _  <- pure vec
--    V.MVector t _  <- pure vec
    (t', Nothing,i) <- lookupRTypeInContext t context
    pure ([t| Ptr $t' |], Nothing,i)

  rev vecConT pt context = do
    AppT vecCon t <- pure pt
    if vecCon /= vecConT
      then mempty
      else do
        t' <- lookupHTypeInContext t context
        pure (Ptr Mutable <$> t' <*> pure ())

pointers :: Q Context
pointers = do
    ptrConT <- [t| Ptr |]
    pure (Context ([rule],[rev ptrConT],[constPtr,mutPtr],"pointers"))
  where
  rule pt context = do
    ff <- pure pt
    case ff of
      FOptional t _ -> do
        (t', Nothing,i) <- lookupRTypeInContext t context
        pure ([t| $t' |], Nothing,i)
      FByReference t _ -> do
        (t', Nothing,i) <- lookupRTypeInContext t context
        pure ([t| Ptr $t' |], Nothing,i)
      _ -> mempty

  rev ptrConT pt context = do
    AppT ptrCon t <- pure pt
    if ptrCon /= ptrConT
      then mempty
      else do
        t' <- lookupHTypeInContext t context
        pure (FByReference <$> t' <*> pure ())

  constPtr = unlines [ "impl<T> MarshalInto<*const T> for *const T {"
                     , "  fn marshal(self) -> *const T { self }"
                     , "}"
                     ]

  mutPtr =   unlines [ "impl<T> MarshalInto<*mut T> for *mut T {"
                     , "  fn marshal(self) -> *mut T { self }"
                     , "}"
                     ]

-- | This maps a Rust function type into the corresponding 'FunPtr' wrapped
-- Haskell function type.
--
-- Note that as a user, you are still responsible for marshalling values of
-- type 'FunPtr'. The reason for this is simple: the GHC runtime has no way of
-- automatically detecting when a pointer to a function is no longer present on
-- the Rust side.
--
-- NOTE: function pointers will not support pointed types that require an intermediate
--       Rust type.

makeIO :: Q HType -> Q HType
makeIO x = do
  ioT <- [t| IO |]
  x1 <- x
  xx <- case x1 of
          AppT i _ -> if i == ioT then x
                                  else [t| IO $x |]
          _ -> [t| IO $x |]
  return xx

--joinPtr :: HType -> HType -> HType
--joinPtr ptrConT a@(AppT p r)
--  | p /= ptrConT = AppT ptrConT a
--  | otherwise = joinPtr ptrConT r
joinPtr :: HType -> HType -> HType
joinPtr ptrConT o@(AppT p2 r@(AppT p3 _))
  | p2 == p3 && ptrConT == p2 = joinPtr ptrConT r
  | otherwise = o
joinPtr _ a = a

functions :: Q Context
functions = do
  funPtrT <- [t| FunPtr |]
  ioT <- [t| IO |]
  ptrConT <- [t| Ptr |]
  pure (Context ([rule ptrConT], [rev funPtrT ioT], [impl],"functions"))
  where
  rule :: HType -> RType -> Context -> First (Q HType, Maybe (Q RType),String)
  rule pt ft context = do
    --BareFn _ C _ (FnDecl args retTy False _) _ <- pure ft
    ff <- pure ft
    case ff of
      (FProcedurePtr _ retTy argTys _) -> do
        args' <- for argTys $ \arg -> do
          (t',_,_) <- lookupRTypeInContext arg context
          pure $ joinPtr pt <$> t'
        retTy' <- do
          (t',_,_) <- lookupRTypeInContext retTy context
          pure t'
        --retTy2 <- pure $ makeIO <$> retTy'
        let retTy2 :: Q HType
            retTy2 = makeIO retTy'
        let hFunTy = foldr (\l r -> [t| $l -> $r |]) retTy2 args'
        let hFunPtr = [t| FunPtr $hFunTy |]
        pure (hFunPtr, Nothing,"FProcedurePtr")
      _ -> mempty

  rev :: Type -> Type -> HType -> Context -> First (Q RType)
  rev funPtrT ioT ft context = do
    AppT funPtr t <- pure ft
    if funPtr /= funPtrT
      then mempty
      else do
        let ts = getApps t
            args = init ts

            ret   = last ts
            ret'  = case ret of
                      AppT io r | io == ioT -> r
                      r -> r
            ret'' = case ret' of
                      TupleT 0 -> Nothing
                      r -> Just r

        argsRs <- traverse (`lookupHTypeInContext` context) args
        retRs <- traverse (`lookupHTypeInContext` context) ret''
        let argsRs' :: Q [Arg ()]
            argsRs' = map (\a -> Arg Nothing a ()) <$> sequence argsRs
--          argsRs2 :: Q [RType]
--          argsRs2 = sequence argsRs
        let decl = FnDecl <$> argsRs' <*> sequence retRs <*> pure False <*> pure ()
        pure (BareFn Normal C [] <$> decl <*> pure ())
--        pure (FProcedurePtr "dummy" retRs argsRs2 <*> pure mempty )
--

  getApps :: Type -> [Type]
  getApps (AppT e1 e2) = e1 : getApps e2
  getApps e = [e]

  -- TODO: this only goes up to 16
  impl :: String
  impl = unlines (map implN [1..16])

    where
    vars = [ l : i | i <- "" : map show [(1 :: Int)..], l <- ['T'..'Z'] ]

    implN :: Int -> String
    implN n = let vs = intercalate "," (take n vars)
                  f = "extern \"C\" fn (" ++ vs ++ ") -> R"
              in unlines [ "impl<" ++ vs ++ ",R> MarshalInto<" ++ f ++ "> for (" ++ f ++ ") {"
                         , "  fn marshal(self) -> (" ++ f ++ ") { self }"
                         , "}"
                         ]

max_stars_repo_name: hasanalrasyid/inline-fortran

max_issues_repo_name: hasanalrasyid/inline-fortran

max_forks_repo_name: hasanalrasyid/inline-fortran


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE StandaloneDeriving #-}

module HLearn.Models.Distributions.Multivariate.MultiNormal
    ( MultiNormal (..)
    )
    where

-- import qualified Control.ConstraintKinds as CK
import qualified Data.Foldable as F
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as VU
import Data.Vector.Unboxed.Deriving

import Control.DeepSeq
import GHC.ST
import GHC.TypeLits

import Foreign.Storable
import Numeric.LinearAlgebra hiding ((<>))
import qualified Numeric.LinearAlgebra as LA

import HLearn.Algebra
import HLearn.Models.Distributions.Common

-------------------------------------------------------------------------------
-- data types

data MultiNormalVec (n::Nat) prob = MultiNormalVec
    { q0 :: !prob
    , q1 :: !(VU.Vector prob)
    , q2 :: !(V.Vector (VU.Vector prob))
    }
    deriving (Read,Show,Eq,Ord)

instance NFData (MultiNormalVec n prob) where
    rnf mn = seq mn ()

newtype MultiNormal prob (xs::[*]) = MultiNormal (MultiNormalVec (Length xs) prob)
    deriving (Read,Show,Eq,Ord,NFData)

deriving instance (Monoid  (MultiNormalVec (Length xs) prob)) => Monoid (MultiNormal prob xs)
deriving instance (Abelian (MultiNormalVec (Length xs) prob)) => Abelian (MultiNormal prob xs)
deriving instance (Group   (MultiNormalVec (Length xs) prob)) => Group (MultiNormal prob xs)
deriving instance (Module  (MultiNormalVec (Length xs) prob)) => Module (MultiNormal prob xs)

-------------------------------------------------------------------------------
-- algebra

instance (Num prob, VU.Unbox prob, SingI n) => Abelian (MultiNormalVec n prob)
instance (Num prob, VU.Unbox prob, SingI n) => Monoid (MultiNormalVec n prob) where
    mempty = MultiNormalVec
        { q0 = 0
        , q1 = VU.replicate n 0
        , q2 = V.replicate n (VU.replicate n 0)
        }
        where
            n = fromIntegral $ fromSing (sing :: Sing n)
    mn1 `mappend` mn2 = MultiNormalVec
        { q0 = (q0 mn1) + (q0 mn2)
        , q1 = VU.zipWith (+) (q1 mn1) (q1 mn2)
        , q2 = V.zipWith (VU.zipWith (+)) (q2 mn1) (q2 mn2)
        }

instance (Num prob, VU.Unbox prob, SingI n) => Group (MultiNormalVec n prob) where
    inverse mn = MultiNormalVec
        { q0 = negate $ q0 mn
        , q1 = VU.map negate (q1 mn)
        , q2 = V.map (VU.map negate) (q2 mn)
        }

instance (Num prob) => HasRing (MultiNormalVec n prob) where
    type Ring (MultiNormalVec n prob) = prob

instance (Num prob, VU.Unbox prob, SingI n) => Module (MultiNormalVec n prob) where
    r .* mn = MultiNormalVec
        { q0 = r * q0 mn
        , q1 = VU.map (r*) (q1 mn)
        , q2 = V.map (VU.map (r*)) (q2 mn)
        }

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

instance (Num prob) => HasRing (MultiNormal prob xs) where
    type Ring (MultiNormal prob xs) = prob

-------------------------------------------------------------------------------
-- training

instance (SingI n, Num prob, VU.Unbox prob) => HomTrainer (MultiNormalVec n prob) where
    type Datapoint (MultiNormalVec n prob) = (VU.Vector prob)
    train1dp dp = MultiNormalVec
        { q0 = 1
        , q1 = dp
        , q2 = V.generate n (\i -> VU.generate n (\j -> (dp VU.! i)*(dp VU.! j)))
        }
        where
            n = fromIntegral $ fromSing (sing :: Sing n)

instance
    ( SingI (Length xs)
    , Num prob
    , VU.Unbox prob
    , HList2List (Datapoint (MultiNormal prob xs)) prob
    ) => HomTrainer (MultiNormal prob xs)
        where
    type Datapoint (MultiNormal prob xs) = HList xs
    train1dp dp = MultiNormal $ train1dp $ VU.fromList $ hlist2list dp

instance (Num prob) => NumDP (MultiNormal prob xs) where
    numdp (MultiNormal mn) = q0 mn

-------------------------------------------------------------------------------
-- distribution

class (Probabilistic dist) => Covariance dist where
    covar :: dist -> Matrix (Probability dist)

instance
    ( VU.Unbox prob
    , SingI k
    , Num prob
    ) => Probabilistic (MultiNormalVec k prob)
        where
    type Probability (MultiNormalVec k prob) = prob

instance
    ( VU.Unbox prob
    , SingI k
    , Fractional prob
    , Enum prob
    , Storable prob
    ) => Covariance (MultiNormalVec k prob)
        where
    covar mn = (k><k) $
            [ (1/(n-1))*( mij - ((q1 mn VU.! j)*((q1 mn) VU.! i))/n )
            | ((i,j),mij) <- assocs
            ]
        where
            assocs = zip [(i,j) | i<-[0..k-1],j<-[0..k-1]] (concat $ V.toList $ fmap VU.toList $ q2 mn)
            mean i = ((q1 mn) VU.! i)/n
            n = q0 mn
            k = fromIntegral $ fromSing (sing :: Sing k)

instance
    ( HList2List (HList dpL) prob
    , VU.Unbox prob
    , Num prob
    , SingI (FromNat1 (Length1 dpL))
    ) => Probabilistic (MultiNormal prob dpL)
        where
    type Probability (MultiNormal prob dpL) = prob

instance
    ( HList2List (HList dpL) prob
    , VU.Unbox prob
    , Floating prob
    , Field prob
    , Enum prob
    , SingI (FromNat1 (Length1 dpL))
--     , Covariance (MultiNormal dpL prob)
    , Storable prob
    ) => PDF (MultiNormal prob dpL)
        where
    pdf (MultiNormal dist) dpL = 1/(sqrt $ (2*pi)^(k)*(det sigma))*(exp $ (-1/2)*(top) )
        where
            top=minElement $ ((trans $ x `sub` mu) LA.<> (inv sigma) LA.<> (x `sub` mu))

            k = fromIntegral $ fromSing (sing :: Sing (Length dpL)) :: Int
            x = k><1 $ hlist2list dpL
            n = q0 dist
            sigma = covar dist
            mu = k><1 $ map (/n) $ VU.toList (q1 $ dist)
--             covarM = (k><k) $ elems $ covar dist

-------------------------------------------------------------------------------
-- tests

-- ds = map (listArray (0,2))
--     [[1,2,4]
--     ,[2,5,6]
--     ,[3,1,1]
--     ]

-- test = train ds :: MultiNormalArray Double 3

ds2 = map VU.fromList
    [[1,2,4]
    ,[2,5,6]
    ,[3,1,1]
    ]

test2 = train ds2 :: MultiNormalVec 3 Double

ds =
    [ 1:::2:::3:::HNil
    , 2:::5:::6:::HNil
    , 3:::1:::1:::HNil
    , 3:::2:::1:::HNil
    ]
test = train ds :: MultiNormal Double '[Double,Double,Double]

max_stars_repo_name: silky/HLearn

max_issues_repo_name: silky/HLearn

max_forks_repo_name: silky/HLearn


text:
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeFamilies #-}
module PCA.Distribution where

import Numeric.LinearAlgebra.Data (Vector, Matrix)

data family Distr m v

data instance Distr Double Double = Gamma Double Double

data instance
     Distr (Vector Double) (Vector Double) = MGamma (Vector Double)
                                                    (Vector Double)

data instance
     Distr (Vector Double) (Matrix Double) = MNormal (Vector Double)
                                                     (Matrix Double)

data instance
     Distr (Matrix Double) (Matrix Double) = MatrixMNormal (Matrix
                                                              Double)
                                                           (Matrix Double)

class Mean m v where
    mean :: Distr m v -> m

instance Mean Double Double where
    mean (Gamma s r) = s / r

instance Mean (Vector Double) (Vector Double) where
    mean (MGamma s r) = s / r

instance Mean (Vector Double) (Matrix Double) where
    mean (MNormal m v) = m

instance Mean (Matrix Double) (Matrix Double) where
    mean (MatrixMNormal m v) = m

class Variance m v  where
    variance :: Distr m v -> v

instance Variance (Vector Double) (Matrix Double) where
    variance (MNormal m v) = v
instance Variance (Matrix Double) (Matrix Double) where
    variance (MatrixMNormal m v) = v

max_stars_repo_name: DbIHbKA/vbpca

max_issues_repo_name: DbIHbKA/vbpca

max_forks_repo_name: DbIHbKA/vbpca