text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE Strict           #-}
module FourierPinwheel.AsteriskGaussian where

import           Data.Array.Repa             as R
import           Data.Complex
import           Data.List                   as L
import           Data.Vector.Generic         as VG
import           Data.Vector.Storable        as VS
import           Data.Vector.Unboxed         as VU
import           DFT.Plan
import           Filter.Utils
import           FourierPinwheel.Array
import           FourierPinwheel.Hypergeo1F1
import           Math.Gamma
import           Pinwheel.FourierSeries2D
import           Utils.Distribution
import           Utils.Parallel

-- The envelope in R2 is r ^ alpha, where -2 < alpha < -0.5
{-# INLINE asteriskGaussian #-}
asteriskGaussian ::
     ( VG.Vector vector (Complex e)
     , RealFloat e
     , NFData (vector (Complex e))
     , Gamma (Complex e)
     , Unbox e
     )
  => Int
  -> Int
  -> e
  -> e
  -> e
  -> e
  -> [vector (Complex e)]
asteriskGaussian numR2Freqs thetaFreq alpha periodR2 periodEnv std =
  parMap
    rdeepseq
    (\angularFreq ->
       let arr =
             R.map
               (* (sqrt (pi / std) *
                   exp (fromIntegral angularFreq ^ 2 / (-4) / std) :+
                   0)) .
             centerHollowArray numR2Freqs $
             analyticalFourierCoefficients1
               numR2Freqs
               1
               angularFreq
               0
               alpha
               periodR2
               periodEnv
        in VG.convert . toUnboxed . computeS $ arr)
    [-thetaFreq .. thetaFreq]

{-# INLINE asteriskGaussianEnvelope #-}
asteriskGaussianEnvelope ::
     (R.Source s (Complex Double))
  => DFTPlan
  -> Int
  -> Int
  -> Double
  -> Double
  -> Double
  -> Double
  -> R.Array s DIM2 (Complex Double)
  -> IO (VS.Vector (Complex Double))
asteriskGaussianEnvelope plan numR2Freqs thetaFreq alpha periodR2 periodEnv stdTheta filter = do
  let planID = DFTPlanID DFT1DG [numR2Freqs, numR2Freqs] [0, 1]
      inversePlanID = DFTPlanID IDFT1DG [numR2Freqs, numR2Freqs] [0, 1]
      vecs =
        asteriskGaussian numR2Freqs thetaFreq alpha periodR2 periodEnv stdTheta
  filterF <-
    dftExecute plan planID . VU.convert . toUnboxed . computeS . makeFilter2D $
    filter
  asteriskGaussianF <- dftExecuteBatchP plan planID vecs
  fmap VS.concat .
    dftExecuteBatchP plan inversePlanID .
    parMap rdeepseq (VS.zipWith (*) filterF) $
    asteriskGaussianF

-- The envelope is r^alpha X Gaussian
asteriskGaussian2 ::
     DFTPlan
  -> Int
  -> Int
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> IO (VS.Vector (Complex Double))
asteriskGaussian2 plan numR2Freqs thetaFreq alpha periodR2 periodEnv stdR2 stdTheta = do
  let centerFreq = div numR2Freqs 2
      gaussian2D =
        fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. i' :. j') ->
          let i = i' - centerFreq
              j = j' - centerFreq
           in exp
                (pi * fromIntegral (i ^ 2 + j ^ 2) /
                 ((-1) * periodR2 ^ 2 * stdR2 ^ 2)) /
              (2 * pi * stdR2 ^ 2) :+
              0
  asteriskGaussianEnvelope
    plan
    numR2Freqs
    thetaFreq
    alpha
    periodR2
    periodEnv
    stdTheta
    gaussian2D

-- The envelope is r^alpha X LowPass
asteriskGaussianLowPass ::
     DFTPlan
  -> Int
  -> Int
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> IO (VS.Vector (Complex Double))
asteriskGaussianLowPass plan numR2Freqs thetaFreq alpha periodR2 periodEnv stdTheta radius = do
  let lpf =
        fromUnboxed (Z :. numR2Freqs :. numR2Freqs) $
        idealLowPassFilter radius periodR2 numR2Freqs
  asteriskGaussianEnvelope
    plan
    numR2Freqs
    thetaFreq
    alpha
    periodR2
    periodEnv
    stdTheta
    lpf

{-# INLINE asteriskGaussianFull #-}
asteriskGaussianFull ::
     ( VG.Vector vector (Complex e)
     , RealFloat e
     , NFData (vector (Complex e))
     , Gamma (Complex e)
     , Unbox e
     )
  => Int
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> e
  -> e
  -> [vector (Complex e)]
asteriskGaussianFull numR2Freqs thetaFreq rFreq alpha periodR2 periodEnv stdTheta stdR =
  let zeroVec = VG.replicate (numR2Freqs ^ 2) 0
   in parMap
        rdeepseq
        (\(radialFreq, angularFreq) ->
           let pinwheel =
                 analyticalFourierCoefficients1
                   numR2Freqs
                   1
                   angularFreq
                   radialFreq
                   alpha
                   periodR2
                   periodEnv
               -- arr =
               --   R.map
               --     (* ((-- gaussian1DFreq (fromIntegral angularFreq) stdTheta *
               --          gaussian1DFourierCoefficients
               --            (fromIntegral radialFreq)
               --            (log periodEnv)
               --            stdR
               --          :+
               --          0) -- *
               --         -- cis
               --         --   (2 * pi / log periodEnv * fromIntegral radialFreq *
               --         --    log 0.5)
               --        )) $
               --   -- if radialFreq == 0 && angularFreq == 0
               --   --   then pinwheel
               --   --   else
                   -- centerHollowArray numR2Freqs pinwheel
               arr =
                 R.map
                   (* ((1 / (1 + (2 * pi / log periodEnv * fromIntegral radialFreq )^2)) :+ 0)
                      ) $
                   centerHollowArray numR2Freqs pinwheel
            in if angularFreq == 0
                 then VG.convert . toUnboxed . computeS $ arr
                 else zeroVec)
        [ (radialFreq, angularFreq)
        | radialFreq <-  [-rFreq .. rFreq]
        , angularFreq <- [-thetaFreq .. thetaFreq]
        ]

{-# INLINE gaussianFull #-}
gaussianFull ::
     ( VG.Vector vector (Complex e)
     , RealFloat e
     , NFData (vector (Complex e))
     , Gamma (Complex e)
     , Unbox e
     )
  => Int
  -> Int
  -> Int
  -> e
  -> e
  -> [vector (Complex e)]
gaussianFull numR2Freqs thetaFreq rFreq periodR2 stdR2 =
  let zeroVec = VG.replicate (numR2Freqs ^ 2) 0
      periodEnv = periodR2 ^ 2 / 4
   in parMap
        rdeepseq
        (\(radialFreq, angularFreq) ->
           if angularFreq == 0 -- && radialFreq == 0
             then VG.convert . toUnboxed . computeS $
                  fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. i' :. j') ->
                    let i = fromIntegral $ i' - div numR2Freqs 2
                        j = fromIntegral $ j' - div numR2Freqs 2
                     in (gaussian2DFourierCoefficients i j periodR2 stdR2 :+ 0) /
                        (1 :+ (2 * pi / log periodEnv * fromIntegral radialFreq)^2)
             else zeroVec)
        [ (radialFreq, angularFreq)
        | radialFreq <- [-rFreq .. rFreq]
        , angularFreq <- [-thetaFreq .. thetaFreq]
        ]

{-# INLINE asteriskGaussianFullEnvelope #-}
asteriskGaussianFullEnvelope ::
     (R.Source s (Complex Double))
  => DFTPlan
  -> Int
  -> Int
  -> Int
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> R.Array s DIM2 (Complex Double)
  -> IO (VS.Vector (Complex Double))
asteriskGaussianFullEnvelope plan numR2Freqs thetaFreq rFreq alpha periodR2 periodEnv stdTheta stdR filter = do
  let planID = DFTPlanID DFT1DG [numR2Freqs, numR2Freqs] [0, 1]
      inversePlanID = DFTPlanID IDFT1DG [numR2Freqs, numR2Freqs] [0, 1]
      vecs =
        asteriskGaussianFull
          numR2Freqs
          thetaFreq
          rFreq
          alpha
          periodR2
          periodEnv
          stdTheta
          stdR
      -- vecs =
      --  gaussianFull
      --      numR2Freqs
      --      thetaFreq
      --      rFreq
      --      periodR2
      --      stdR
  filterF <-
    dftExecute plan planID . VU.convert . toUnboxed . computeS . makeFilter2D $
    filter
  asteriskGaussianF <- dftExecuteBatchP plan planID vecs
  fmap VS.concat .
    dftExecuteBatchP plan inversePlanID .
    parMap rdeepseq (VS.zipWith (*) filterF) $
    asteriskGaussianF
  -- return . VS.concat $ vecs

-- The envelope is r^alpha X LowPass
asteriskGaussianFullLowPass ::
     DFTPlan
  -> Int
  -> Int
  -> Int
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> IO (VS.Vector (Complex Double))
asteriskGaussianFullLowPass plan numR2Freqs thetaFreq rFreq alpha periodR2 periodEnv stdTheta stdR radius = do
  let lpf =
        fromUnboxed (Z :. numR2Freqs :. numR2Freqs) $
        idealLowPassFilter radius periodR2 numR2Freqs
  asteriskGaussianFullEnvelope
    plan
    numR2Freqs
    thetaFreq
    rFreq
    alpha
    periodR2
    periodEnv
    stdTheta
    stdR
    lpf

-- The envelope is r^alpha X LowPass
asteriskGaussian2Full ::
     DFTPlan
  -> Int
  -> Int
  -> Int
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> IO (VS.Vector (Complex Double))
asteriskGaussian2Full plan numR2Freqs thetaFreq rFreq alpha periodR2 periodEnv stdTheta stdR stdR2 = do
  let centerFreq = div numR2Freqs 2
      gaussian2D =
        fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. i' :. j') ->
          let i = fromIntegral $ i' - centerFreq
              j = fromIntegral $ j' - centerFreq
           in gaussian2DFourierCoefficients i j periodR2 stdR2 :+ 0
  asteriskGaussianFullEnvelope
    plan
    numR2Freqs
    thetaFreq
    rFreq
    alpha
    periodR2
    periodEnv
    stdTheta
    stdR
    gaussian2D

asteriskGaussianRFull ::
     DFTPlan
  -> Int
  -> Int
  -> Int
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> IO (VS.Vector (Complex Double))
asteriskGaussianRFull plan numR2Freqs thetaFreq rFreq alpha periodR2 periodEnv stdTheta stdR stdR2 = do
  let planID = DFTPlanID DFT1DG [numR2Freqs, numR2Freqs] [0, 1]
      inversePlanID = DFTPlanID IDFT1DG [numR2Freqs, numR2Freqs] [0, 1]
      centerFreq = div numR2Freqs 2
      gaussian2D =
        fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. i' :. j') ->
          let i = i' - centerFreq
              j = j' - centerFreq
           in exp
                (pi * fromIntegral (i ^ 2 + j ^ 2) /
                 ((-1) * periodR2 ^ 2 * stdR2 ^ 2)) /
              (2 * pi * stdR2 ^ 2) :+
              0
      r =
        centerHollowArray' numR2Freqs $ analyticalFourierCoefficients3 numR2Freqs 1 0 0 alpha periodR2 periodEnv
      vecs =
        asteriskGaussianFull
          numR2Freqs
          thetaFreq
          rFreq
          alpha
          periodR2
          periodEnv
          stdTheta
          stdR
  gaussian2DF <-
    dftExecute plan planID . VU.convert . toUnboxed . computeS . makeFilter2D $
    gaussian2D
  rF <-
    dftExecute plan planID . VU.convert . toUnboxed . computeS . makeFilter2D $
    r
  let filterF = VS.zipWith (\x y -> x + y ^ 2) gaussian2DF rF
  asteriskGaussianF <- dftExecuteBatchP plan planID vecs
  fmap VS.concat .
    dftExecuteBatchP plan inversePlanID .
    parMap rdeepseq (VS.zipWith (*) filterF) $
    asteriskGaussianF


text:
{-# LANGUAGE Arrows #-}

-- Homogeneous (or depth-preserving) functions over perfectly balanced trees.

module Main where

import Control.Arrow
import Control.Category
import Data.Complex
import Prelude hiding (id, (.))

infixr 4 :&:

-- Consider the following non-regular type of perfectly balanced trees,
-- or `powertrees' (cf Jayadev Misra's powerlists):

data Pow a = Zero a | Succ (Pow (Pair a))
	deriving Show

type Pair a = (a, a)

-- Here are some example elements:

tree0 = Zero 1
tree1 = Succ (Zero (1, 2))
tree2 = Succ (Succ (Zero ((1, 2), (3, 4))))
tree3 = Succ (Succ (Succ (Zero (((1, 2), (3, 4)), ((5, 6), (7, 8))))))

-- The elements of this type have a string of constructors expressing
-- a depth n as a Peano numeral, enclosing a nested pair tree of 2^n
-- elements.  The type definition ensures that all elements of this type
-- are perfectly balanced binary trees of this form.  (Such things arise
-- in circuit design, eg Ruby, and descriptions of parallel algorithms.)
-- And the type system will ensure that all legal programs preserve
-- this structural invariant.
--
-- The only problem is that the type constraint is too restrictive, rejecting
-- many of the standard operations on these trees.  Typically you want to
-- split a tree into two subtrees, do some processing on the subtrees and
-- combine the results.  But the type system cannot discover that the two
-- results are of the same depth (and thus combinable).  We need a type
-- that says a function preserves depth.  Here it is:

data Hom a b = (a -> b) :&: Hom (Pair a) (Pair b)

-- A homogeneous (or depth-preserving) function is an infinite sequence of
-- functions of type Pair^n a -> Pair^n b, one for each depth n.  We can
-- apply a homogeneous function to a powertree by selecting the function
-- for the required depth:

apply :: Hom a b -> Pow a -> Pow b
apply (f :&: fs) (Zero x) = Zero (f x)
apply (f :&: fs) (Succ t) = Succ (apply fs t)

-- Having defined apply, we can forget about powertrees and do all our
-- programming with Hom's.  Firstly, Hom is an arrow:

instance Category Hom where
	id = id :&: id
	(f :&: fs) . (g :&: gs) = (f . g) :&: (fs . gs)

instance Arrow Hom where
	arr f = f :&: arr (f *** f)
	first (f :&: fs) =
		first f :&: (arr transpose >>> first fs >>> arr transpose)

transpose :: ((a,b), (c,d)) -> ((a,c), (b,d))
transpose ((a,b), (c,d)) = ((a,c), (b,d))

-- arr maps f over the leaves of a powertree.

-- The composition >>> composes sequences of functions pairwise.
--
-- The *** operator unriffles a powertree of pairs into a pair of powertrees,
-- applies the appropriate function to each and riffles the results.
-- It defines a categorical product for this arrow category.

-- When describing algorithms, one often provides a pure function for the
-- base case (trees of one element) and a (usually recursive) expression
-- for trees of pairs.

-- For example, a common divide-and-conquer pattern is the butterfly, where
-- one recursive call processes the odd-numbered elements and the other
-- processes the even ones (cf Geraint Jones and Mary Sheeran's Ruby papers):

butterfly :: (Pair a -> Pair a) -> Hom a a
butterfly f = id :&: proc (x, y) -> do
		x' <- butterfly f -< x
		y' <- butterfly f -< y
		returnA -< f (x', y')

-- The recursive calls operate on halves of the original tree, so the
-- recursion is well-defined.

-- Some examples of butterflies:

rev :: Hom a a
rev = butterfly swap
	where	swap (x, y) = (y, x)

unriffle :: Hom (Pair a) (Pair a)
unriffle = butterfly transpose

-- Batcher's sorter for bitonic sequences:

bisort :: Ord a => Hom a a
bisort = butterfly cmp
	where	cmp (x, y) = (min x y, max x y)

-- This can be used (with rev) as the merge phase of a merge sort.
--
sort :: Ord a => Hom a a
sort = id :&: proc (x, y) -> do
		x' <- sort -< x
		y' <- sort -< y
		yr <- rev -< y'
		p <- unriffle -< (x', yr)
		bisort2 -< p
	where _ :&: bisort2 = bisort

-- Here is the scan operation, using the algorithm of Ladner and Fischer:

scan :: (a -> a -> a) -> a -> Hom a a
scan op b = id :&: proc (x, y) -> do
		y' <- scan op b -< op x y
		l <- rsh b -< y'
		returnA -< (op l x, y')

-- The auxiliary function rsh b shifts each element in the tree one place to
-- the right, placing b in the now-vacant leftmost position, and discarding
-- the old rightmost element:

rsh :: a -> Hom a a
rsh b = const b :&: proc (x, y) -> do
		w <- rsh b -< y
		returnA -< (w, x)

-- Finally, here is the Fast Fourier Transform:

type C = Complex Double

fft :: Hom C C
fft = id :&: proc (x, y) -> do
		x' <- fft -< x
		y' <- fft -< y
		r <- roots (-1) -< ()
		let z = r*y'
		unriffle -< (x' + z, x' - z)

-- The auxiliary function roots r (where r is typically a root of unity)
-- populates a tree of size n (necessarily a power of 2) with the values
-- 1, w, w^2, ..., w^(n-1), where w^n = r.

roots :: C -> Hom () C
roots r = const 1 :&: proc _ -> do
		x <- roots r' -< ()
		unriffle -< (x, x*r')
	where	r' = if imagPart s >= 0 then -s else s
		s = sqrt r

-- Miscellaneous functions:

rrot :: Hom a a
rrot = id :&: proc (x, y) -> do
		w <- rrot -< y
		returnA -< (w, x)

ilv :: Hom a a -> Hom (Pair a) (Pair a)
ilv f = proc (x, y) -> do
		x' <- f -< x
		y' <- f -< y
		returnA -< (x', y')

scan' :: (a -> a -> a) -> a -> Hom a a
scan' op b = proc x -> do
		l <- rsh b -< x
		(id :&: ilv (scan' op b)) -< op l x

riffle :: Hom (Pair a) (Pair a)
riffle = id :&: proc ((x1, y1), (x2, y2)) -> do
		x <- riffle -< (x1, x2)
		y <- riffle -< (y1, y2)
		returnA -< (x, y)

invert :: Hom a a
invert = id :&: proc (x, y) -> do
		x' <- invert -< x
		y' <- invert -< y
		unriffle -< (x', y')

carryLookaheadAdder :: Hom (Bool, Bool) Bool
carryLookaheadAdder = proc (x, y) -> do
		carryOut <- rsh (Just False) -<
			if x == y then Just x else Nothing
		Just carryIn <- scan plusMaybe Nothing -< carryOut
		returnA -< x `xor` y `xor` carryIn
	where	plusMaybe x Nothing = x
		plusMaybe x (Just y) = Just y
		False `xor` b = b
		True `xor` b = not b

-- Global conditional for SIMD

ifAll :: Hom a b -> Hom a b -> Hom (a, Bool) b
ifAll fs gs = ifAllAux snd (arr fst >>> fs) (arr fst >>> gs)
	where	ifAllAux :: (a -> Bool) -> Hom a b -> Hom a b -> Hom a b
		ifAllAux p (f :&: fs) (g :&: gs) =
			liftIf p f g :&: ifAllAux (liftAnd p) fs gs
		liftIf p f g x = if p x then f x else g x
		liftAnd p (x, y) = p x && p y

maybeAll :: Hom a c -> Hom (a, b) c -> Hom (a, Maybe b) c
maybeAll (n :&: ns) (j :&: js) =
	choose :&: (arr dist >>> maybeAll ns (arr transpose >>> js))
	where	choose (a, Nothing) = n a
		choose (a, Just b) = j (a, b)
		dist ((a1, b1), (a2, b2)) = ((a1, a2), zipMaybe b1 b2)
		zipMaybe (Just x) (Just y) = Just (x, y)
		zipMaybe _ _ = Nothing

main = do
	print (apply rev tree3)
	print (apply invert tree3)
	print (apply (invert >>> sort) tree3)
	print (apply (scan (+) 0) tree3)

ext: hs

max_stars_repo_name: ocharles/ghc

max_issues_repo_name: ocharles/ghc

max_forks_repo_name: ocharles/ghc


text:
-- | Common infrastructure for CVode/ARKode
{-# OPTIONS_GHC -Wno-name-shadowing #-}
module Numeric.Sundials.Common where

import Foreign.C.Types
import Foreign.Ptr
import Foreign.Storable (peek, poke)
import Foreign.C.String
import Numeric.Sundials.Types
import qualified Numeric.Sundials.Foreign as T
import qualified Data.Vector as V
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VSM
import Data.Maybe
import Numeric.LinearAlgebra.HMatrix hiding (Vector)
import GHC.Prim
import Control.Monad.IO.Class
import Control.Monad.Cont
import Control.Exception
import Katip
import Data.Aeson
import qualified Data.Text as T
import qualified Data.Text.Encoding as T
import qualified Data.ByteString as BS
import Control.Monad.Reader
import GHC.Generics (Generic)
import Foreign.ForeignPtr

-- | A collection of variables that we allocate on the Haskell side and
-- pass into the C code to be filled.
data CVars vec = CVars
  { c_diagnostics :: vec SunIndexType
    -- ^ Mutable vector to which we write diagnostic data while
    -- solving. Its size corresponds to the number of fields in
    -- 'SundialsDiagnostics'.
  , c_root_info :: vec CInt
    -- ^ Just a temporary vector (of the size equal to the number of event
    -- specs) that we use to get root info. Isn't used for output.
  , c_event_index :: vec CInt
    -- ^ For each event occurrence, this indicates which of the events
    -- occurred. Size: max_num_events.
  , c_event_time :: vec CDouble
    -- ^ For each event occurrence, this indicates the time of the
    -- occurrence. Size: max_num_events.
  , c_n_events :: vec CInt
    -- ^ Vector of size 1 that gives the total number of events occurred.
  , c_n_rows :: vec CInt
    -- ^ The total number of rows in the output matrix.
  , c_output_mat :: vec CDouble
    -- ^ The output matrix stored in the row-major order.
    -- Dimensions: (1 + dim) * (2 * max_events + nTs).
  , c_actual_event_direction :: vec CInt
    -- ^ Vector of size max_num_events that gives the direction of the
    -- occurred event.
  , c_local_error :: vec CDouble
    -- ^ Vector containing local error estimates. Size: the dimensionality
    -- of the system.
  , c_var_weight :: vec CDouble
    -- ^ Vector containing variable weights (derived from the tolerances).
    -- Size: the dimensionality of the system.
  , c_local_error_set :: vec CInt
    -- The flag (size 1) indicating whether c_local_error is filled with meaningful
    -- values. *Should be initialized with 0.*
  }

allocateCVars :: OdeProblem -> IO (CVars (VS.MVector RealWorld))
allocateCVars OdeProblem{..} = do
  let dim = VS.length odeInitCond
  c_diagnostics <- VSM.new 11
  c_root_info <- VSM.new $ V.length odeEvents
  c_event_index <- VSM.new odeMaxEvents
  c_event_time <- VSM.new odeMaxEvents
  c_actual_event_direction <- VSM.new odeMaxEvents
  c_n_events <- VSM.new 1
  c_n_rows <- VSM.new 1
  c_local_error <- VSM.new dim
  c_var_weight <- VSM.new dim
  c_local_error_set <- VSM.new 1
  VSM.write c_local_error_set 0 0
  c_output_mat <- VSM.new $
    (1 + dim) * (2 * odeMaxEvents + VS.length odeSolTimes)
  return CVars {..}

-- NB: the mutable CVars must not be used after this
freezeCVars :: CVars (VS.MVector RealWorld) -> IO (CVars VS.Vector)
freezeCVars CVars{..} = do
  c_diagnostics <- VS.unsafeFreeze c_diagnostics
  c_root_info <- VS.unsafeFreeze c_root_info
  c_event_index <- VS.unsafeFreeze c_event_index
  c_event_time <- VS.unsafeFreeze c_event_time
  c_actual_event_direction <- VS.unsafeFreeze c_actual_event_direction
  c_n_events <- VS.unsafeFreeze c_n_events
  c_n_rows <- VS.unsafeFreeze c_n_rows
  c_output_mat <- VS.unsafeFreeze c_output_mat
  c_local_error <- VS.unsafeFreeze c_local_error
  c_var_weight <- VS.unsafeFreeze c_var_weight
  c_local_error_set <- VS.unsafeFreeze c_local_error_set
  return CVars {..}

-- | Similar to 'CVars', except these are immutable values that are
-- accessed (read-only) by the C code and specify the system to be solved.
data CConsts = CConsts
  { c_dim :: SunIndexType -- ^ the dimensionality (number of variables/equations)
  , c_method :: CInt -- ^ the ODE method (specific to the solver)
  , c_n_sol_times :: CInt
  , c_sol_time :: VS.Vector CDouble
  , c_init_cond :: VS.Vector CDouble
  , c_rhs :: FunPtr OdeRhsCType
  , c_rhs_userdata :: Ptr UserData
  , c_rtol :: CDouble
  , c_atol :: VS.Vector CDouble
  , c_n_event_specs :: CInt
  , c_event_fn :: CDouble -> Ptr T.SunVector -> Ptr CDouble -> Ptr () -> IO CInt
  , c_apply_event
      :: CInt -- number of triggered events
      -> Ptr CInt -- event indices
      -> CDouble -- time
      -> Ptr T.SunVector -- y
      -> Ptr T.SunVector -- new y
      -> Ptr CInt -- (out) stop the solver?
      -> Ptr CInt -- (out) record the event?
      -> IO CInt
  , c_jac_set :: CInt
  , c_jac :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix
          -> Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector
          -> IO CInt
  , c_sparse_jac :: CInt
      -- ^ If 0, use a dense matrix.
      -- If non-0, use a sparse matrix with that number of non-zero
      -- elements.
  , c_requested_event_direction :: VS.Vector CInt
  , c_next_time_event :: IO CDouble
  , c_max_events :: CInt
  , c_minstep :: CDouble
  , c_fixedstep :: CDouble
  , c_max_n_steps :: SunIndexType
  , c_max_err_test_fails :: CInt
  , c_init_step_size_set :: CInt
  , c_init_step_size :: CDouble
  }

data Solver = CVode | ARKode
  deriving Show

class Method method where
  methodToInt :: method -> CInt
  methodSolver :: Solver

withCConsts
  :: Method method
  => ODEOpts method
  -> OdeProblem
  -> (CConsts -> IO r)
  -> IO r
withCConsts ODEOpts{..} OdeProblem{..} = runContT $ do
  let
    dim = VS.length c_init_cond
    c_init_cond = coerce odeInitCond
    c_dim = fromIntegral dim
    c_n_sol_times = fromIntegral . VS.length $ odeSolTimes
    c_sol_time = coerce odeSolTimes
    c_rtol = relTolerance odeTolerances
    c_atol = either (VS.replicate dim) id $ absTolerances odeTolerances
    c_minstep = coerce minStep
    c_fixedstep = coerce fixedStep
    c_max_n_steps = fromIntegral maxNumSteps
    c_max_err_test_fails = fromIntegral maxFail
    c_init_step_size_set = fromIntegral . fromEnum $ isJust initStep
    c_init_step_size = coerce . fromMaybe 0 $ initStep
    c_n_event_specs = fromIntegral $ V.length odeEvents
    c_requested_event_direction = V.convert $ V.map (directionToInt . eventDirection) odeEvents
    c_event_fn t y_ptr out_ptr _ptr = do
      y <- sunVecVals <$> peek y_ptr
      let vals = V.convert $ V.map (\ev -> coerce (eventCondition ev) t y) odeEvents
      -- FIXME: We should be able to use poke somehow
      T.vectorToC vals (fromIntegral c_n_event_specs) out_ptr
      return 0
    c_apply_event n_events event_indices_ptr t y_ptr y'_ptr stop_solver_ptr record_event_ptr = do
      event_indices <- vecFromPtr event_indices_ptr (fromIntegral n_events)
      -- Apparently there's no safe version of
      y_vec <- peek y_ptr
      EventHandlerResult{..} <-
        odeEventHandler
          (coerce t :: Double)
          (coerce $ sunVecVals y_vec :: VS.Vector Double)
          (VS.map fromIntegral event_indices :: VS.Vector Int)
      poke y'_ptr $ SunVector
        { sunVecN = sunVecN y_vec
        , sunVecVals = coerce eventNewState
        }
      poke stop_solver_ptr . fromIntegral $ fromEnum eventStopSolver
      poke record_event_ptr . fromIntegral $ fromEnum eventRecord
      return 0
    c_max_events = fromIntegral odeMaxEvents
    c_next_time_event = coerce odeTimeBasedEvents
    c_jac_set = fromIntegral . fromEnum $ isJust odeJacobian
    c_jac t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
      case odeJacobian of
        Nothing   -> undefined
        Just jacI -> do j <- matrixToSunMatrix . jacI (coerce t) <$> (coerce $ sunVecVals <$> peek y)
                        case jacobianRepr of
                          DenseJacobian -> poke jacS j
                          SparseJacobian spat -> poke (castPtr jacS) (T.SparseMatrix spat j)
                        return 0
    c_sparse_jac = case jacobianRepr of
      SparseJacobian (T.SparsePattern spat) -> VS.sum (VS.map fromIntegral spat)
      DenseJacobian -> 0
    c_method = methodToInt odeMethod

  (c_rhs, c_rhs_userdata) <-
    case odeRhs of
      OdeRhsC ptr u -> return (ptr, u)
      OdeRhsHaskell fun -> do
        let
          funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr UserData -> IO CInt
          funIO t y f _ptr = do
            sv <- peek y
            r <- fun t (sunVecVals sv)
            poke f $ SunVector { sunVecN = sunVecN sv
                               , sunVecVals = r
                               }
            return 0
        funptr <- ContT $ bracket (mkOdeRhsC funIO) freeHaskellFunPtr
        return (funptr, nullPtr)
  return CConsts{..}

matrixToSunMatrix :: Matrix Double -> T.SunMatrix
matrixToSunMatrix m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs }
  where
    nr = fromIntegral $ rows m
    nc = fromIntegral $ cols m
    vs = coerce . VS.concat $ toColumns m

-- Contrary to the documentation, it appears that CVodeGetRootInfo
-- may use both 1 and -1 to indicate a root, depending on the
-- direction of the sign change. See near the end of cvRootfind.
intToDirection :: Integral d => d -> Maybe CrossingDirection
intToDirection d =
  case d of
    1  -> Just Upwards
    -1 -> Just Downwards
    _  -> Nothing

-- | Almost inverse of 'intToDirection'. Map 'Upwards' to 1, 'Downwards' to
-- -1, and 'AnyDirection' to 0.
directionToInt :: Integral d => CrossingDirection -> d
directionToInt d =
  case d of
    Upwards -> 1
    Downwards -> -1
    AnyDirection -> 0

foreign import ccall "wrapper"
  mkOdeRhsC :: OdeRhsCType -> IO (FunPtr OdeRhsCType)

assembleSolverResult
  :: OdeProblem
  -> CInt
  -> CVars VS.Vector
  -> IO (Either ErrorDiagnostics SundialsSolution)
assembleSolverResult OdeProblem{..} ret CVars{..} = do
  let
    dim = VS.length odeInitCond
    n_rows = fromIntegral . VS.head $ c_n_rows
    output_mat = coerce . reshape (dim + 1) . subVector 0 ((dim + 1) * n_rows) $ c_output_mat
    (local_errors, var_weights) =
      if c_local_error_set VS.! 0 == 0
        then (mempty, mempty)
        else coerce (c_local_error, c_var_weight)
    diagnostics = SundialsDiagnostics
      (fromIntegral $ c_diagnostics VS.!0)
      (fromIntegral $ c_diagnostics VS.!1)
      (fromIntegral $ c_diagnostics VS.!2)
      (fromIntegral $ c_diagnostics VS.!3)
      (fromIntegral $ c_diagnostics VS.!4)
      (fromIntegral $ c_diagnostics VS.!5)
      (fromIntegral $ c_diagnostics VS.!6)
      (fromIntegral $ c_diagnostics VS.!7)
      (fromIntegral $ c_diagnostics VS.!8)
      (fromIntegral $ c_diagnostics VS.!9)
      (toEnum . fromIntegral $ c_diagnostics VS.! 10)
  return $
    if ret == T.cV_SUCCESS
      then
        Right $ SundialsSolution
          { actualTimeGrid = extractTimeGrid output_mat
          , solutionMatrix = dropTimeGrid output_mat
          , diagnostics
          }
      else
        Left ErrorDiagnostics
          { partialResults = output_mat
          , errorCode = fromIntegral ret
          , errorEstimates = local_errors
          , varWeights = var_weights
          }
  where
    -- The time grid is the first column of the result matrix
    extractTimeGrid :: Matrix Double -> VS.Vector Double
    extractTimeGrid = head . toColumns
    dropTimeGrid :: Matrix Double -> Matrix Double
    dropTimeGrid = fromColumns . tail . toColumns

-- | The common solving logic between ARKode and CVode
solveCommon
  :: (Method method, Katip m)
  => (CConsts -> CVars (VS.MVector RealWorld) -> ReportErrorFn -> IO CInt)
      -- ^ the CVode/ARKode solving function; mostly inline-C code
  -> ODEOpts method
  -> OdeProblem
  -> m (Either ErrorDiagnostics SundialsSolution)
solveCommon solve_c opts problem@(OdeProblem{..})

  | VS.null odeInitCond = -- 0-dimensional (empty) system

    return . Right $ SundialsSolution
      { actualTimeGrid = odeSolTimes
      , solutionMatrix = (VS.length odeSolTimes >< 0) []
      , diagnostics = emptyDiagnostics
      }

  | otherwise = do

    report_error <- logWithKatip
    liftIO $ do -- the rest is in the IO monad
    vars <- allocateCVars problem
    ret <- withCConsts opts problem $ \consts ->
      solve_c consts vars report_error
    frozenVars <- freezeCVars vars
    assembleSolverResult problem ret frozenVars

-- | An auxiliary function to construct a storable vector from a C pointer
-- and length.
--
-- There doesn't seem to be a safe version of 'VS.unsafeFromForeignPtr0',
-- nor a way to clone an immutable vector, so we emulate it via an
-- intermediate mutable vector.
vecFromPtr
  :: VS.Storable a
  => Ptr a
  -> Int
  -> IO (VS.Vector a)
vecFromPtr ptr n = do
  fptr <- newForeignPtr_ ptr
  let mv = VSM.unsafeFromForeignPtr0 fptr n
  VS.freeze mv -- this does the copying and makes the whole thing safe

----------------------------------------------------------------------
--                           Logging
----------------------------------------------------------------------

-- | The Katip payload for logging Sundials errors
data SundialsErrorContext = SundialsErrorContext
  { sundialsErrorCode :: !Int
  , sundialsErrorModule :: !T.Text
  , sundialsErrorFunction :: !T.Text
  } deriving Generic
instance ToJSON SundialsErrorContext
instance ToObject SundialsErrorContext
instance LogItem SundialsErrorContext where
  payloadKeys _ _ = AllKeys

type ReportErrorFn =
  (  CInt    -- error code
  -> CString -- module name
  -> CString -- function name
  -> CString -- the message
  -> Ptr ()  -- user data (ignored)
  -> IO ()
  )

logWithKatip
  :: Katip m
  => m ReportErrorFn
logWithKatip = do
  log_env <- getLogEnv
  return $
    \err_code c_mod_name c_func_name c_msg _userdata ->
    -- See Note [CV_TOO_CLOSE]
    if err_code == T.cV_TOO_CLOSE then pure () else do
    let
      toText :: CString -> IO T.Text
      toText = fmap T.decodeUtf8 . BS.packCString
    mod_name <- toText c_mod_name
    func_name <- toText c_func_name
    msg <- toText c_msg
    let
      severity :: Severity
      severity =
        if err_code <= 0
          then ErrorS
          else WarningS
      errCtx :: SundialsErrorContext
      errCtx = SundialsErrorContext
        { sundialsErrorCode = fromIntegral err_code
        , sundialsErrorModule = mod_name
        , sundialsErrorFunction = func_name
        }
    flip runReaderT log_env . unKatipT $ do
      logF errCtx "sundials" severity (logStr msg)

max_stars_repo_name: tweag/hmatrix-sundials

max_issues_repo_name: tweag/hmatrix-sundials

max_forks_repo_name: haskell-numerics/hmatrix-sundials


text:
module Numeric.FFT.Execute ( execute ) where

import           Control.Monad               (when)
import qualified Control.Monad               as CM
import           Control.Monad.ST
import           Data.Complex
import qualified Data.IntMap.Strict          as IM
import           Data.STRef
import qualified Data.Vector                 as V
import           Data.Vector.Unboxed
import qualified Data.Vector.Unboxed.Mutable as MV
import           Prelude                     hiding (concatMap, foldr, length,
                                              map, mapM_, null, reverse, sum,
                                              zip, zipWith)

import           Numeric.FFT.Special
import           Numeric.FFT.Types
import           Numeric.FFT.Utils

-- | Main FFT plan execution driver.
execute :: Plan -> Direction -> VCD -> VCD
execute (Plan dlinfo perm base) dir h
  | n == 1 = h
  | V.null dlinfo = runST $ do
                          mhin <- case perm of
                            Nothing -> thaw h
                            Just p  -> unsafeThaw $ backpermute h p
                          mhout <- MV.replicate n 0
                          applyBase base sign mhin mhout
                          when (dir == Inverse) $ do
                            let s = 1.0 / fromIntegral n :+ 0
                            CM.forM_ [0..n-1] $ \i -> do
                              x <- MV.unsafeRead mhout i
                              MV.unsafeWrite mhout i $ s * x
                          unsafeFreeze mhout
  | otherwise=  fullfft
  where
    n = length h              -- Input vector length.
    bsize = baseSize base     -- Size of base transform.

    -- Root of unity sign.
    sign = case dir of
      Forward -> 1
      Inverse -> -1

    -- Apply Danielson-Lanczos steps and base transform to digit
    -- reversal ordered input vector.
    fullfft = runST $ do
      mhin <- case perm of
            Nothing -> thaw h
            Just p  -> unsafeThaw $ backpermute h p
      mhtmp <- MV.replicate n 0
      multBase mhin mhtmp
      mhr <- newSTRef (mhtmp, mhin)
      V.forM_ dlinfo $ \dlstep -> do
        (mh0, mh1) <- readSTRef mhr
        dl sign dlstep mh0 mh1
        writeSTRef mhr (mh1, mh0)
      mhs <- readSTRef mhr
      let vout = fst mhs
      when (dir == Inverse) $ do
        let s = 1.0 / fromIntegral n :+ 0
        CM.forM_ [0..n-1] $ \i -> do
          x <- MV.unsafeRead vout i
          MV.unsafeWrite vout i $ s * x
      unsafeFreeze vout

    -- Multiple base transform application for "bottom" of algorithm.
    multBase :: MVCD s -> MVCD s -> ST s ()
    multBase xmin xmout =
      V.zipWithM_ (applyBase base sign)
                  (slicemvecs bsize xmin) (slicemvecs bsize xmout)

-- | Monadic FFT plan execution driver -- used by Rader's algorithm
-- for convolutions.
executeM :: Plan -> Direction -> MVCD s -> MVCD s -> ST s ()
executeM (Plan dlinfo perm base) dir hin hout =
  if n == 1
  then MV.copy hout hin
  else do
    htmp <- MV.replicate n 0

    -- Input permutation.
    case perm of
      Nothing -> MV.copy htmp hin
      Just p  -> backpermuteM n p hin htmp

    -- Apply Danielson-Lanczos steps and base transform to digit
    -- reversal ordered input vector.
    multBase htmp hout
    mhr <- newSTRef (hout, htmp)
    V.forM_ dlinfo $ \dlstep -> do
      (mh0, mh1) <- readSTRef mhr
      dl sign dlstep mh0 mh1
      writeSTRef mhr (mh1, mh0)
    when (odd $ V.length dlinfo) $ MV.copy hout htmp

    -- Output scaling for inverse transform.
    when (dir == Inverse) $ do
      let s = 1.0 / fromIntegral n :+ 0
      forM_ (enumFromN 0 n) $ \i -> do
        x <- MV.unsafeRead hout i
        MV.unsafeWrite hout i $ s * x
  where
    n = MV.length hin         -- Input vector length.
    bsize = baseSize base     -- Size of base transform.

    -- Root of unity sign.
    sign = case dir of
      Forward -> 1
      Inverse -> -1

    -- Multiple base transform application for "bottom" of algorithm.
    multBase :: MVCD s -> MVCD s -> ST s ()
    multBase xmin xmout =
      V.zipWithM_ (applyBase base sign)
                  (slicemvecs bsize xmin) (slicemvecs bsize xmout)

-- | Single Danielson-Lanczos step: process all duplicates and
-- concatenate into a single vector.
dl :: Int -> (Int, Int, VVVCD, VVVCD) -> MVCD s -> MVCD s -> ST s ()
dl sign (wfac, split, dmatp, dmatm) mhin mhout =
  V.zipWithM_ doone (slicemvecs wfac mhin) (slicemvecs wfac mhout)
  where
    -- Twiddled diagonal entries in row r, column c (both
    -- zero-indexed), where each row and column if a wfac x wfac
    -- matrix.
    dmat = if sign == 1 then dmatp else dmatm
    d r c = (dmat V.! r) V.! c

    -- Size of each diagonal sub-matrix.
    ns = wfac `div` split

    -- Index vectors.
    nsidxs = enumFromN 0 ns
    splitidxs = enumFromN 1 (split-1)

    -- Process one duplicate by processing all rows and writing the
    -- results into a single output vector.
    doone :: MVCD s -> MVCD s -> ST s ()
    doone vin vout = do
      let vs = (slicemvecs ns vin, slicemvecs ns vout)
      mapM_ (single vs) $ enumFromN 0 split
      where
        -- Multiply a single block by its appropriate diagonal
        -- elements and accumulate the result.
        mult :: VMVCD s -> MVCD s -> Int -> Bool -> Int -> ST s ()
        mult vins vo r first c = do
          let vi = vins V.! c
              dvals = d r c
          forM_ nsidxs $ \i -> do
            xi <- MV.unsafeRead vi i
            xo <- if first then return 0 else MV.unsafeRead vo i
            MV.unsafeWrite vo i (xo + xi * dvals ! i)
        -- Multiply all blocks by the corresponding diagonal
        -- elements in a single row.
        single :: (VMVCD s, VMVCD s) -> Int -> ST s ()
        single (vis, vos) r = do
          mult vis (vos V.! r) r True 0
          mapM_ (mult vis (vos V.! r) r False) splitidxs
        -- single (vis, vos) r =
        --   let m = mult vis (vos V.! r) r
        --   in do
        --     m True 0
        --     mapM_ (m False) splitidxs

-- | Apply a base transform to a single vector.
applyBase :: BaseTransform -> Int -> MVCD s -> MVCD s -> ST s ()

-- Simple DFT algorithm.
applyBase (DFTBase sz wsfwd wsinv) sign mhin mhout = do
  h <- freeze mhin
  forM_ (enumFromN 0 sz) $ \i -> MV.unsafeWrite mhout i (doone h i)
  where ws = if sign == 1 then wsfwd else wsinv
        doone h i = sum $ zipWith (*) h $
                    generate sz (\k -> ws ! (i * k `mod` sz))

-- Special hard-coded cases.
applyBase (SpecialBase sz) sign mhin mhout =
  case IM.lookup sz specialBases of
    Just f  -> f sign mhin mhout
    Nothing -> error "invalid problem size for SpecialBase"

-- Rader prime-length FFT.
applyBase (RaderBase sz outperm bfwd binv csz cplan) sign mhin mhout = do
  -- Padding size.
  let pad = csz - (sz - 1)

  -- Permuted input vector padded to next greater power of two size
  -- for fast convolution.
  apad <- MV.replicate csz 0
  forM_ (enumFromN 0 csz) $ \i -> do
    val <- if i == 0 then MV.unsafeRead mhin 1
           else if i > pad
                then MV.unsafeRead mhin $ i - pad + 1
                else return 0
    MV.unsafeWrite apad i val

  -- FFT-based convolution calculation.
  convtmp <- MV.replicate csz 0
  executeM cplan Forward apad convtmp
  let bmult = if sign == 1 then bfwd else binv
  forM_ (enumFromN 0 csz) $ \i -> do
    x <- MV.unsafeRead convtmp i
    MV.unsafeWrite convtmp i $ x * (bmult ! i)
  executeM cplan Inverse convtmp apad
  conv <- unsafeFreeze apad

  -- Input vector sum.
  sumhref <- newSTRef 0
  forM_ (enumFromN 0 sz) $ \i -> do
    val <- MV.unsafeRead mhin i
    modifySTRef sumhref (+ val)
  sumh <- readSTRef sumhref

  -- Write output based on output generator index ordering.
  h0 <- MV.unsafeRead mhin 0
  forM_ (enumFromN 0 sz) $ \i -> do
    let (idx, val) = case i of
          0 -> (0, sumh)
          _ -> (outperm ! (i - 1), h0 + conv ! (i - 1))
    MV.unsafeWrite mhout idx val

max_stars_repo_name: ian-ross/arb-fft

max_issues_repo_name: ian-ross/arb-fft

max_forks_repo_name: ian-ross/arb-fft


text:
{-# LANGUAGE GADTs, OverloadedStrings, DeriveFunctor #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Data.Matrix.MatrixMarket.Internal
-- Copyright   :  (c) Marco Zocca 2017
-- License     :  BSD2 (see the file LICENSE)
--
-- Maintainer  :  zocca marco gmail
-- Stability   :  experimental
-- Portability :  portable
--
-- Attoparsec parser and serializer for the NIST MatrixMarket format. The parser logic originally appeared in `accelerate-examples` and it is reused here (courtesy of T.McDonell and the `accelerate` developers) with some amendments.
--
-- In this version:
-- * Numbers are represented with Scientific notation instead of floating point
-- * Parsing rules are a bit relaxed to accommodate various whitespace corner cases
--
-----------------------------------------------------------------------------
module Data.Matrix.MatrixMarket.Internal
       (readMatrix, readMatrix', readArray,
        writeMatrix, writeMatrix', writeArray,
        Matrix(..), Array(..),
        Format (Coordinate, Array), Structure (General, Symmetric, Hermitian, Skew),
        nnz, dim, numDat,
        dimArr, numDatArr,
        ImportError) where

import Control.Applicative                      hiding ( many )

import Data.Int
import qualified Data.Char as C
import Data.Complex
import qualified Data.Scientific as S
import Data.Attoparsec.ByteString.Char8 hiding (I)
-- import qualified Data.ByteString as BS
import qualified Data.ByteString.Lazy.Char8 as B
import qualified Data.Attoparsec.Lazy           as L
import qualified Data.ByteString.Lazy           as L

import Control.Monad.Catch
import Control.Exception.Common

-- | Specifies either sparse or dense storage.  In sparse (\"coordinate\")
-- storage, elements are given in (i,j,x) triplets for matrices (or (i,x) for
-- vectors).  Indices are 1-based, so that A(1,1) is the first element of a
-- matrix, and x(1) is the first element of a vector.
--
-- In dense (\"array\") storage, elements are given in column-major order.
--
-- In both cases, each element is given on a separate line.
--
data Format = Coordinate | Array
    deriving (Eq, Show)

-- | Specifies the element type.  Pattern matrices do not have any elements,
-- only indices, and only make sense for coordinate matrices and vectors.
--
data Field  = R | C | I | P
    deriving (Eq, Show)

-- | Specifies any special structure in the matrix.  For symmetric and hermitian
-- matrices, only the lower-triangular part of the matrix is given. For skew
-- matrices, only the entries below the diagonal are stored.
--
data Structure = General | Symmetric | Hermitian | Skew
    deriving (Eq, Show)

-- | Sparse matrix in coordinate form (row, column, entry)
-- NB: indices are 1-based i.e. A(1,1) is the top-left entry of matrix A
data Matrix a = RMatrix (Int, Int) Int Structure [(Int, Int, a)]
              | CMatrix (Int, Int) Int Structure [(Int, Int, Complex a)]
              | PatternMatrix (Int,Int) Int  Structure [(Int32,Int32)]
              | IntMatrix (Int,Int) Int Structure [(Int32,Int32,Int)]
              deriving (Eq, Show)

-- | Array, i.e. a DENSE matrix (also used to represent vectors as n-by-1 matrices)
data Array a = RArray (Int, Int) Structure [a]
             | CArray (Int, Int) Structure [Complex a]
             deriving (Eq, Show)

--------------------------------------------------------------------------------
-- Attoparsec combinators
--------------------------------------------------------------------------------

comment :: Parser ()
comment = char '%' *> skipWhile (not . eol) *> endOfLine
  where
    eol w = w `elem` ("\n\r" :: String)

floating :: Parser S.Scientific
floating = skipSpace' *> scientific

integral :: Integral a => Parser a
integral = skipSpace' *> decimal

format :: Parser Format
format =  string "coordinate" *> pure Coordinate
      <|> string "array"      *> pure Array
      <?> "matrix format"

field :: Parser Field
field =  string "real"    *> pure R
     <|> string "complex" *> pure C
     <|> string "integer" *> pure I
     <|> string "pattern" *> pure P
     <?> "matrix field"

structure :: Parser Structure
structure =  string "general"        *> pure General
         <|> string "symmetric"      *> pure Symmetric
         <|> string "hermitian"      *> pure Hermitian
         <|> string "skew-symmetric" *> pure Skew
         <?> "matrix structure"

header :: Parser (Format,Field,Structure)
header =  string "%%MatrixMarket matrix"
       >> (,,) <$> (skipSpace' *> format)
               <*> (skipSpace' *> field)
               <*> (skipSpace' *> structure)
               <*  endOfLine
               <?> "MatrixMarket header"

extentMatrix :: Parser (Int,Int,Int)
extentMatrix = do
  [m,n,l] <- skipSpace' *> count 3 integral <* endOfLine
  return (m,n,l)

extentArray :: Parser (Int,Int)
extentArray = do
  [m,n] <- skipSpace' *> count 2 integral <* endOfLine
  return (m,n)

line3 :: Integral i => Parser a -> Parser (i,i,a)
line3 f = (,,) <$> integral
               <*> integral
               <*> f
               <*  (endOfLine <|> L.endOfInput)

skipSpace' :: Parser String
skipSpace' = many' space

--------------------------------------------------------------------------------
-- Load and parse
--------------------------------------------------------------------------------

matrix :: Parser (Matrix S.Scientific)
matrix = do
  (f, t, s) <- header
  (m, n, l) <- skipMany comment *> extentMatrix
  if f /= Coordinate
    then fail "matrix is not in Coordinate format"
    else
    case t of
     R -> RMatrix (m,n) l s <$> many1 (line3 floating)
     C -> CMatrix (m,n) l s <$> many1 (line3 ((:+) <$> floating <*> floating))
     I -> IntMatrix     (m,n) l s <$> many1 (line3 integral)
     P -> PatternMatrix (m,n) l s <$> many1 ((,) <$> integral <*> integral)

array :: Parser (Array S.Scientific)
array = do
  (f, t, s) <- header
  (m, n) <- skipMany comment *> extentArray
  if f /= Array
    then fail "array is not in Array format"
    else
     case t of
       R -> RArray (m,n) s <$> many1 floating
       C -> CArray (m,n) s <$> many1 ((:+) <$> floating <*> floating)
       _ -> fail "integer and pattern cases not relevant for the dense case"

-- | Load a matrix (sparse, i.e. in Coordinate format) from file
readMatrix :: FilePath -> IO (Matrix S.Scientific)
readMatrix file = L.readFile file >>= readMatrix'

-- | Load a matrix (sparse, i.e. in Coordinate format) from a bytestring.
readMatrix' :: L.ByteString -> IO (Matrix S.Scientific)
readMatrix' chunks =
  case L.parse matrix chunks of
    L.Fail _ _ msg      -> throwM (FileParseError "readMatrix" msg)
    L.Done _ mtx        -> return mtx

-- | Load a dense matrix (i.e. a matrix or vector in Array format) from file
readArray :: FilePath -> IO (Array S.Scientific)
readArray file = do
  chunks <- L.readFile file
  case L.parse array chunks of
    L.Fail _ _ msg      -> throwM (FileParseError "readArray" msg)
    L.Done _ mtx        -> return mtx

--------------------------------------------------------------------------------
-- Write to file
--------------------------------------------------------------------------------

showFormat :: Format -> String
showFormat = map C.toLower <$> show

showField :: Field -> String
showField f = case f of R -> "real"
                        C -> "complex"
                        I -> "integer"
                        P -> "pattern"

showStruct :: Structure -> String
showStruct = map C.toLower <$> show

-- %%MatrixMarket matrix coordinate real general

headerStr :: Format -> Field -> Structure -> L.ByteString
headerStr f t s =
    B.pack $ unwords ["%%MatrixMarket matrix",
                      showFormat f, showField t, showStruct s]

nl :: L.ByteString
nl = toLBS "\n"

showLines :: (a -> String) -> [a] -> L.ByteString
showLines showf d = L.concat (L.pack . withNewline . showf <$> d) where
  withNewline x = toEnum . C.ord <$> x ++ "\n"

-- | Serialize a sparse matrix in Coordinate format
writeMatrix :: Show b => FilePath -> Matrix b -> IO ()
writeMatrix file = L.writeFile file . writeMatrix'

-- | Serialize a sparse matrix in Coordinate format as a bytestring
writeMatrix' :: Show b => Matrix b -> L.ByteString
writeMatrix' mat =
  case mat of (RMatrix d nz s dat) ->
                matrixByteString d nz R s dat
              (CMatrix d nz s dat) ->
                matrixByteString d nz C s dat
              (IntMatrix d nz s dat) ->
                matrixByteString d nz I s dat
              _ -> error "writeMatrix : PatternMatrix not implemented yet"
     where
       matrixByteString di nz t s d =
         L.concat [headerStr Coordinate t s,
                   nl,
                   headerSzMatrix di nz,
                   nl,
                   showLines sf3 d]
         where
         sf3 (i,j,x) = unwords [show i, show j, show x]
         headerSzMatrix (m,n) numz = B.pack $ unwords [show m, show n, show numz]

-- | Serialize a dense matrix in Array format
writeArray :: Show a => FilePath -> Array a -> IO ()
writeArray file arr =
  case arr of (RArray d s dat) -> L.writeFile file (arrayByteString d R s dat)
              (CArray d s dat) -> L.writeFile file (arrayByteString d C s dat)
    where
      arrayByteString di t s d =
        L.concat [headerStr Array t s,
            nl,
            headerSzArray di,
            nl,
            showLines show d] where
        headerSzArray (m,n) = B.pack $ unwords [show m, show n]

-- | helpers

-- | Number of matrix nonzeros
nnz :: Matrix t -> Int
nnz m = case m of (RMatrix _ nz _ _) -> nz
                  (CMatrix _ nz _ _) -> nz
                  (PatternMatrix _ nz _ _) -> nz
                  (IntMatrix _ nz _ _) -> nz

-- | Matrix size : number of rows, number of columns
dim :: Matrix t -> (Int, Int)
dim m = case m of (RMatrix d _ _ _) -> d
                  (CMatrix d _ _ _) -> d
                  (PatternMatrix d _ _ _) -> d
                  (IntMatrix d _ _ _) -> d

-- | Length of data vector internal to the Matrix; this is _not_ necessarily the actual number of matrix entries because symmetric entries are not stored
numDat :: Matrix t -> Int
numDat m = case m of (RMatrix _ _ _ d) -> length d
                     (CMatrix _ _ _ d) -> length d
                     (PatternMatrix _ _ _ d) -> length d
                     (IntMatrix _ _ _ d) -> length d

-- | Array size : number of rows, number of columns
dimArr :: Array t -> (Int, Int)
dimArr a = case a of (RArray d _ _) -> d
                     (CArray d _ _) -> d

-- | Length of data vector internal to the Array; this is _not_ necessarily the actual number of matrix entries because symmetric entries are not stored
numDatArr :: Array a -> Int
numDatArr a = case a of (RArray _ _ ll) -> length ll
                        (CArray _ _ ll) -> length ll

-- | String -> ByteString
-- NB: do not abuse
-- toBS :: String -> BS.ByteString
-- toBS x = BS.pack $ (toEnum . C.ord) <$> x

-- | String -> lazy ByteString
-- NB: do not abuse
toLBS :: String -> L.ByteString
toLBS x = L.pack $ (toEnum . C.ord) <$> x

-- t1 = toBS "%%MatrixMarket matrix coordinate real general\n5  5  8\n1     1   1.000e+00\n2     2   1.050e+01\n3     3   1.500e-02\n1     4   6.000e+00\n4     2   2.505e+02\n4     4  -2.800e+02\n4     5   3.332e+01\n5     5   1.200e+01\n"

-- t2 = toBS "%%MatrixMarket matrix coordinate real general\n5  5  8\n1     1   1.000e+00\n"


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE AllowAmbiguousTypes #-}

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

{-|
Module      : Numeric.LinearAlgebra.Static.RList
Description : List of static vectors (`R`s)
Copyright   : (c) Ryan Orendorff, 2020
License     : BSD3
Stability   : experimental
-}
module Numeric.LinearAlgebra.Static.RList
    ( RList(..)
    , LList(..)
    , Sum
    , exRList
    , concat
    , blockDiag
    , mapRListSquare
    , applyLList
    )
where

import           GHC.TypeLits
import           Numeric.LinearAlgebra.Static
import Prelude hiding (concat)

-- TODO: Make a fold function, then this can be abstracted.
-- | Sum a type level list of naturals
type family Sum (xs :: [Nat]) where
    Sum '[] = 0
    Sum (n ': ns) = n + Sum ns

-- | A non-empty list of `R` vectors of different sizes.
-- The list is non-empty since we cannot create a vector of size 0 in a
-- reasonable way when using Data.Vector as the backend (as HMatrix).
data RList (ns :: [Nat]) where
    ROne :: (KnownNat n) => R n -> RList '[n]
    (:::) :: (KnownNat n) => R n -> RList ns -> RList (n ': ns)

infixr 6 :::

-- | An example of how to construct an RList.
exRList :: RList '[2, 3]
exRList = konst 1 ::: ROne (konst 2)

-- There is no null vector. Also the HMatrix library does the wrong thing
-- with this example.
--
-- λ> (konst 4 :: R 2) # (konst 0 :: R 0)
-- (vector [4.0,4.0,0.0] :: R 2)

-- | Concatenate an RList of vectors
--
-- Thank you Christiaan Baaij for figuring out how to do concat! This is
-- basically his solution to the problem, with a modification to work with
-- the non-empty version of RList. See the following gist for more details:
-- https://gist.github.com/ryanorendorff/d05c378b71829e3c0c33de462cb9a973
concat :: RList ns -> R (Sum ns)
concat rs = case concatWorker rs of
  RWorker rsW -> rsW

data RWorker n where
  RWorker :: KnownNat n => R n -> RWorker n

concatWorker :: RList ns -> RWorker (Sum ns)
concatWorker (ROne r) = RWorker r
concatWorker (r ::: rs) = case concatWorker rs of
  RWorker rsW -> RWorker (r # rsW)

-- | A non-empty list of `L` matrices of different sizes.
-- The list is non-empty since we cannot create a matrix of size 0 in a
-- reasonable way when using Data.Vector as the backend (as HMatrix).
data LList (ms :: [Nat]) (ns :: [Nat]) where
    LOne :: (KnownNat m, KnownNat n) => L m n -> LList '[m] '[n]
    (:-:) :: (KnownNat m, KnownNat n) => L m n -> LList ms ns
                                      -> LList (m ': ms) (n ': ns)

-- | Concatenate a LList of matrices as a block diagonal.
blockDiag :: LList ms ns -> L (Sum ms) (Sum ns)
blockDiag ls = case blockDiagWorker ls of
  LWorker lsW -> lsW

data LWorker m n where
  LWorker :: (KnownNat m, KnownNat n) => L m n -> LWorker m n

blockDiagWorker :: LList ms ns -> LWorker (Sum ms) (Sum ns)
blockDiagWorker (LOne l) = LWorker l
blockDiagWorker (l :-: ls) = case blockDiagWorker ls of
  LWorker lsW -> LWorker (l ||| konst 0 === konst 0 ||| lsW)

-- | Map some function that preserves the dimension of the vector (a square
-- matrix) over a list of vectors.
mapRListSquare :: (forall n. KnownNat n => R n -> R n)
               -> RList ns -> RList ns
mapRListSquare f (ROne r) = ROne (f r)
mapRListSquare f (r ::: rs) = f r ::: mapRListSquare f rs

-- I have not proved that the LList and RList are the same length, so GHC
-- complains that I have non-exhaustive patterns that cannot be constructed.
-- TODO: Add comment on how this would work in Agda to show how this type
-- of error would be avoided.
-- | Apply a list of static matrices over a list of static vectors
applyLList :: LList ms ns -> RList ns -> RList ms
applyLList (LOne _) (_ ::: _) = error "This case is impossible"
applyLList (_ :-: _) (ROne _) = error "This case is impossible"
applyLList (LOne l) (ROne r) = ROne (l #> r)
applyLList (l :-: ls) (r ::: rs) = (l #> r) ::: applyLList ls rs

-- This function is more complicated, may need to reflect out the position
-- we are in the list in order to write this function.
-- splitRList :: R (Sum ns) -> RList ns
-- splitRList vec = let (r, rs) = split vec in r ::: (splitRList rs)

max_stars_repo_name: ryanorendorff/convex

max_issues_repo_name: ryanorendorff/convex

max_forks_repo_name: ryanorendorff/convex


text:
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Matrix.Eigen
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Eigenvalue decompositions.
--
module Numeric.LinearAlgebra.Matrix.Eigen (
    hermEigen,
    hermEigenvalues,
    hermEigenM_,
    hermEigenvaluesM_,
    ) where

import Control.Monad( when )
import Control.Monad.ST( ST, runST, unsafeIOToST )
import Foreign( allocaArray, nullPtr )
import Text.Printf( printf )

import Foreign.LAPACK( LAPACK, EigJob(..), EigRange(..) )
import qualified Foreign.LAPACK as LAPACK

import Numeric.LinearAlgebra.Types( Herm(..) )
import Numeric.LinearAlgebra.Matrix.Base( Matrix )
import Numeric.LinearAlgebra.Matrix.STBase( STMatrix )
import qualified Numeric.LinearAlgebra.Matrix.STBase as M
import Numeric.LinearAlgebra.Vector( Vector, STVector )
import qualified Numeric.LinearAlgebra.Vector as V

-- | Compute the eigenvalues and eigenvectors of a Hermitian matrix.
-- Return the eigenvalues are in ascending order in the result vector;
-- store the corresponding eigenvectors are in the columns of the result
-- matrix.
hermEigen :: (LAPACK e) => Herm Matrix e -> (Vector Double, Matrix e)
hermEigen (Herm uplo a) = runST $ do
    (_,n) <- M.getDim a
    ma <- M.newCopy a
    mw <- V.new_ n
    mz <- M.new_ (n,n)
    hermEigenM_ (Herm uplo ma) mw mz
    w <- V.unsafeFreeze mw
    z <- M.unsafeFreeze mz
    return (w,z)

-- | Return the eigenvalues of a Hermitian matrix in ascending order.
hermEigenvalues :: (LAPACK e) => Herm Matrix e -> Vector Double
hermEigenvalues (Herm uplo a) = V.create $ do
    (_,n) <- M.getDim a
    ma <- M.newCopy a
    w <- V.new_ n
    hermEigenvaluesM_ (Herm uplo ma) w
    return w

-- | Compute and copy the eigenvalues and eigenvectors of a mutable
-- Hermitian matrix.  This destroys the original Hermitian matrix.
hermEigenM_ :: (LAPACK e)
                  => Herm (STMatrix s) e
                  -> STVector s Double
                  -> STMatrix s e
                  -> ST s ()
hermEigenM_ (Herm uplo a) w z = do
    (ma,na) <- M.getDim a
    (mz,nz) <- M.getDim z
    nw <- V.getDim w
    let n = na

    when (ma /= na) $ error $
        printf ("hermEigenM_"
                ++ " (Herm _ <matrix with dim (%d,%d)>): nonsquare matrix"
               ) ma na
    when ((not . and) [ (ma,na) == (n,n)
                      , nw == n
                      , (mz,nz) == (n,n)
                      ]) $ error $
        printf ("hermEigenM_"
                ++ " (Herm _ <matrix with dim (%d,%d)>)"
                ++ " <vector with dim %d>:"
                ++ " <matrix with dim (%d,%d)>"
                ++ " dimension mismatch")
               ma na nw mz nz

    unsafeIOToST $
        M.unsafeWith a $ \pa lda ->
        V.unsafeWith w $ \pw ->
        M.unsafeWith z $ \pz ldz ->
        allocaArray (2*n) $ \psuppz -> do
            _m <- LAPACK.heevr jobz range uplo n pa lda abstol pw pz ldz
                               psuppz
            return ()
  where
    jobz = EigVectors
    range = AllEigs
    abstol = 0

-- | Compute and copy the eigenvalues of a mutable Hermitian matrix.  This
-- destroys the original Hermitian matrix.
hermEigenvaluesM_ :: (LAPACK e)
                        => Herm (STMatrix s) e
                        -> STVector s Double
                        -> ST s ()
hermEigenvaluesM_ (Herm uplo a) w = do
    (ma,na) <- M.getDim a
    nw <- V.getDim w
    let n = na

    when (ma /= na) $  error $
        printf ("hermEigenvaluesM_"
                ++ " (Herm _ <matrix with dim (%d,%d)>): nonsquare matrix"
               ) ma na

    when ((not . and) [ (ma,na) == (n,n)
                      , nw == n
                      ]) $ error $
        printf ("hermEigenvaluesM_"
                ++ " (Herm _ <matrix with dim (%d,%d)>)"
                ++ " <vector with dim %d>: dimension mismatch")
               ma na nw

    unsafeIOToST $
        M.unsafeWith a $ \pa lda ->
        V.unsafeWith w $ \pw -> do
            _m <- LAPACK.heevr jobz range uplo n pa lda abstol pw pz ldz
                               psuppz
            return ()

  where
    jobz = NoEigVectors
    range = AllEigs
    abstol = 0
    pz = nullPtr
    ldz = 1
    psuppz = nullPtr

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-|
Module: MachineLearning.Random
Description: Random generation utility functions.
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX

Randon generation uitility functions.
-}

module MachineLearning.Random
(
  sample
  , sampleM
  , getRandomRListM
  , getRandomRVectorM
  , getRandomRMatrixM
)

where

import Control.Monad (when)
import System.Random (RandomGen, Random)
import qualified Control.Monad.ST as ST
import qualified Data.Vector as V
import qualified Data.Vector.Storable as SV
import Numeric.LinearAlgebra ((><))
import qualified Data.Vector.Mutable as MV
import qualified Control.Monad.Random as RndM
import MachineLearning.Types (R, Vector, Matrix)

-- | Samples `n` (given as a second parameter) values from `list` (given as a third parameter).
sample :: RandomGen g => g -> Int -> V.Vector a -> (V.Vector a, g)
sample gen n xs = RndM.runRand (sampleM n xs) gen

-- | Samples `n` (given as a second parameter) values from `list` (given as a third parameter) inside RandomMonad.
sampleM :: RandomGen g => Int -> V.Vector a -> RndM.Rand g (V.Vector a)
sampleM n xs = do  -- Random Monad starts
  let rangeList = V.fromList $ zip (repeat 0) [n..(length xs)-1]
  rnds <- randomsInRangesM rangeList
  let (pre, post) = V.splitAt n xs
  let ys = ST.runST $ do  -- ST Monad starts
        mv <- V.thaw pre
        V.zipWithM_ (\val r -> when (r < n) $ MV.write mv (mod r n) val) post rnds
        V.unsafeFreeze mv
  return ys

-- | Returns a list of random values distributed in a closed interval `range`
getRandomRListM :: (RandomGen g, Random a) =>
                   Int             -- ^ list's lengths
                   -> (a, a)       -- ^ range
                   -> RndM.Rand g [a]          -- ^ list of random values inside RandomMonad
getRandomRListM size range = mapM (\_ -> RndM.getRandomR range) [1..size]

-- | Returns a vector of random values distributed in a closed interval `range`
getRandomRVectorM :: RandomGen g =>
                    Int                     -- ^ vector's length
                    -> (R, R)               -- ^ range
                    -> RndM.Rand g Vector   -- vector of randon values inside RandomMonad
getRandomRVectorM size range = SV.fromList <$> getRandomRListM size range

-- | Returns a matrix of random values distributed in a closed interval `range`
getRandomRMatrixM :: RandomGen g =>
                    Int                     -- ^ number of rows
                    -> Int                  -- ^ number of columns
                    -> (R, R)               -- ^ range
                    -> RndM.Rand g Matrix   -- vector of randon values inside RandomMonad
getRandomRMatrixM r c range = (r><c) <$> getRandomRListM (r*c) range

-- | Takes a list of ranges `(lo, hi)`,
-- returns a list of random values uniformly distributed in the list of closed intervals [(lo, hi)].
randomsInRangesM :: (RndM.RandomGen g, RndM.Random a) => V.Vector (a, a) -> RndM.Rand g (V.Vector a)
randomsInRangesM rangeList = mapM RndM.getRandomR rangeList

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Core.Language where

import Control.Monad.Except
import Data.Array
import Data.Complex
import Data.IORef
import Data.Ratio
import System.IO
import Text.Printf (printf)

import Core.Error

-- SCHEME DATA TYPES -----------------------------------------------------------

type SchemeEnv = IORef [(String, IORef SchemeVal)]

type SchemeFun = [SchemeVal] -> ThrowsError SchemeVal

data SchemeVal = Atom String
               | Boolean Bool
               | Char Char
               | String String
               | Integer Int
               | Rational Rational
               | Real Double
               | Complex (Complex Double)
               | List [SchemeVal]
               | ImproperList [SchemeVal] SchemeVal
               | Vector (Array Int SchemeVal)
               | Bytevector (Array Int Integer)
               | Primitive SchemeFun
               | Func { args :: [String],
                        vararg :: Maybe String,
                        body :: [SchemeVal],
                        closure :: SchemeEnv }
               | IOPrimitive ([SchemeVal] -> IOThrowsError SchemeVal)
               | IOPort Handle
               | Unspecified

-- STRING REPRESENTATIONS ------------------------------------------------------

instance Show SchemeVal where show = showVal

showVal (Atom name) = name
showVal (Boolean bool) = if bool then "#t" else "#f"
showVal (Char ch) = show ch
showVal (String str) = printf "\"%s\"" str
showVal (Integer n) = show n
showVal (Rational q) = printf "%d/%d" (numerator q) (denominator q)
showVal (Real x) = show x
showVal (Complex z) = printf "%f+%fi" (realPart z) (imagPart z)
showVal (List xs) = printf "(%s)" (unwordsList xs)
showVal (ImproperList xs x) = printf "(%s . %s)" (unwordsList xs) (showVal x)
showVal (Vector array) = printf "#(%s)" (unwordsArray array)
showVal (Bytevector array) = printf "#u8(%s)" (unwordsArray array)
showVal (Primitive _) = "<primitive>"
showVal (Func args vararg _ _) = showLambda args vararg
showVal (IOPrimitive _) = "<IO primitive>"
showVal (IOPort _) = "<IO port>"
showVal Unspecified = "<unspecified>"

-- UTILS -----------------------------------------------------------------------

-- throws a wrong number of arguments error
throwWrongNumArgs :: Int -> [SchemeVal] -> ThrowsError SchemeVal
throwWrongNumArgs = (throwError .) . (. unwordsList) . WrongNumArgs

-- throws a type mismatch error
throwTypeMismatch :: String -> SchemeVal -> ThrowsError SchemeVal
throwTypeMismatch = (throwError .) . (. show) . TypeMismatch

unwordsList :: Show a => [a] -> String
unwordsList = unwords . map show

unwordsArray :: Show a => Array Int a -> String
unwordsArray = unwordsList . elems

showLambda :: [String] -> Maybe String -> String
showLambda = (. maybe [] (" . " ++)) . printf "(lambda (%s%s) ...)" . unwords

max_stars_repo_name: calintat/scheme

max_issues_repo_name: calintat/scheme

max_forks_repo_name: calintat/scheme


text:
{-# LANGUAGE BangPatterns #-}

import Data.Char
import Data.List
import Debug.Trace
import Numeric.LinearAlgebra hiding (Matrix, Vector)
import qualified Numeric.LinearAlgebra as LA
import System.Random.Mersenne

type FpType = Double
type Vector = LA.Vector FpType
type Matrix = LA.Matrix FpType

-- Initializer

initVector :: MTGen -> Int -> IO (Vector)
initVector rng n = do
    rs <- randoms rng
    let v = fromList (take n rs)
    return v

initMatrix :: MTGen -> Int -> Int -> IO (Matrix)
initMatrix rng n m = do
    rs <- randoms rng :: IO [FpType]
    let v = (n><m) $ take (n*m) rs
    return v

initConnections :: MTGen -> (Layer, Layer) -> IO (WeightMatrix)
initConnections rng (l₁,l₂) = do
    w <- initMatrix rng (neurons l₂) (neurons l₁)
    return w

initMlp :: MTGen -> [Layer] -> IO (Network)
initMlp rng ls = do
    ws <- mapM (initConnections rng) $ zip ls (tail ls)
    let net = mlp ls ws
    return net

-- Helper

vecToStr :: Vector -> String
vecToStr v = "[ " ++ intercalate "  " (map show (toList v)) ++ " ]"

vecsToStr :: [Vector] -> String
vecsToStr vs = intercalate "\n" (map vecToStr vs)

datasetToStr :: Dataset -> String
datasetToStr d = "dataset: " ++ show l ++ " elements, " ++ show u ++ " targets\n"
    where u = nub (map snd d)
          l = length d

-- Layers

--TODO: allow 2d tensors for activation as well to facilitate softmax
data Layer = Layer {
    ϕ         :: FpType -> FpType,
    ϕ'        :: FpType -> FpType,
    neurons   :: Int,
    layerType :: String
}

instance Show Layer where
    show (Layer { neurons = n, layerType = t }) = t ++ ", " ++ show n ++ " neurons"

linearLayer :: Int -> Layer
linearLayer n = Layer {
    ϕ         = id,
    ϕ'        = id,
    neurons   = n,
    layerType = "linear"
}

hyperbolicLayer :: Int -> Layer
hyperbolicLayer n = Layer {
    ϕ         = tanh,
    ϕ'        = \x -> 1 - tanh² x,
    neurons   = n,
    layerType = "tanh"
} where tanh² = tanh . tanh

sigmoidLayer :: Int -> Layer
sigmoidLayer n = Layer {
    ϕ         = sigmoid,
    ϕ'        = \x -> sigmoid x * (1 - sigmoid x),
    neurons   = n,
    layerType = "sigmoid"
} where sigmoid x = 0.5 * (1 + tanh (0.5 * x))

-- Networks

data Network = Network {
    layers  :: [Layer],
    weights :: [WeightMatrix],
    netType :: String
}

type WeightMatrix = Matrix

--TODO: It would be nicer to enforce #layers == #weights+1 by zipping.
mlp :: [Layer] -> [WeightMatrix] -> Network
mlp ls ws = debug Network { layers = ls, weights = ws, netType = "mlp" }
    where text  = map (\x -> "Layer " ++ show (fst x) ++ ": " ++ show (snd x)) (zip [1..] ls)
          debug = trace $ "Creating MLP\n" ++ intercalate "\n" text ++ "\n"

forward :: Vector -> (Layer, WeightMatrix) -> Vector
forward x (l,w) = mapVector (ϕ l) $ w <> x

activate :: Network -> Vector -> [Vector]
activate (Network { layers = ls, weights = ws }) x = debug result
    where result = scanl forward ϕx $ zip (tail ls) ws
          ϕx     = mapVector ϕ₁ x
          ϕ₁     = ϕ (head ls)
          debug  = trace ("activate input\n" ++ vecToStr x ++ "\n\nactivate output\n" ++ vecsToStr result ++ "\n")

-- Datasets

--TODO: decouple classification/regression datasets properly, add bias s.t. it works flawlessly in all usages
type Dataset = [(Vector, Vector)]

-- Training

data Trainer = Trainer {
    η         :: FpType,
    terminate :: TrainState -> Bool
}

data TrainState = TrainState {
    nEpochs :: Int,
    δs      :: [FpType]
} deriving (Show)

type Error = FpType

--TODO: more useful termination criteria
nEpochTrainer :: FpType -> Int -> Trainer
nEpochTrainer lr n = debug Trainer { η = lr, terminate = (> n) . nEpochs }
    where debug = trace ("Creating nEpochTrainer\n" ++ "η = " ++ show lr ++ "\nmaxEpochs = " ++ show n ++ "\n")

--TODO: actually perform epoch inside ST monad, put actual result
epoch :: Dataset -> Network -> Trainer -> (Network, Error)
epoch d net t = (net, 0)

_train :: TrainState -> Dataset -> Network -> Trainer -> Network
_train s d net t = do
    let (net',δ') = epoch d net t
    let s' = update s where update s = TrainState { nEpochs = 1 + (nEpochs s), δs = δs s ++ [δ'] }
    let statusMessage = "δ_" ++ show (nEpochs s') ++ " = " ++ show δ' ++ "\n" ++ show s' ++ "\n"

    if (terminate t) s'
        then trace ("Finished Training.\n") net'
        else trace statusMessage $ _train s' d net' t

--TODO: support testing set/cross validation
train :: Dataset -> Network -> Trainer -> Network
train = _train s₀ where s₀ = TrainState { nEpochs = 0, δs = [] }

-- Evaluation

--TODO

-- Testing

main :: IO()
main = do
    rng <- newMTGen Nothing

    let l₁ = linearLayer 5
    let l₂ = sigmoidLayer 3
    let l₃ = linearLayer 5
    let layers = [l₁,l₂,l₃]

    !net <- initMlp rng layers
    xs <- mapM (initVector rng) (take 100 $ repeat 5)

    let ys = xs
    let dataset = zip xs ys
    let trainer = nEpochTrainer 0.1 5
    let !trainedNet = train dataset net trainer

    putStrLn $ "Done!\n"


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE DefaultSignatures #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE Trustworthy #-}
---------------------------------------------------------------------------
-- |
-- Copyright   :  (C) 2012-2015 Edward Kmett
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  experimental
-- Portability :  non-portable
--
-- Simple matrix operation for low-dimensional primitives.
---------------------------------------------------------------------------
module Linear.Trace
  ( Trace(..)
  , frobenius
  ) where

import Control.Monad as Monad
import Linear.V0
import Linear.V1
import Linear.V2
import Linear.V3
import Linear.V4
import Linear.Plucker
import Linear.Quaternion
import Linear.V
import Linear.Vector
import Data.Complex
import Data.Distributive
import Data.Foldable as Foldable
import Data.Functor.Bind as Bind
import Data.Functor.Compose
import Data.Functor.Product
import Data.Hashable
import Data.HashMap.Lazy
import Data.IntMap (IntMap)
import Data.Map (Map)

-- $setup
-- >>> import Data.Complex
-- >>> import Debug.SimpleReflect.Vars
-- >>> import Linear.V2

class Functor m => Trace m where
  -- | Compute the trace of a matrix
  --
  -- >>> trace (V2 (V2 a b) (V2 c d))
  -- a + d
  trace :: Num a => m (m a) -> a
#ifndef HLINT
  default trace :: (Foldable m, Num a) => m (m a) -> a
  trace = Foldable.sum . diagonal
  {-# INLINE trace #-}
#endif

  -- | Compute the diagonal of a matrix
  --
  -- >>> diagonal (V2 (V2 a b) (V2 c d))
  -- V2 a d
  diagonal :: m (m a) -> m a
#ifndef HLINT
  default diagonal :: Monad m => m (m a) -> m a
  diagonal = Monad.join
  {-# INLINE diagonal #-}
#endif

instance Trace IntMap where
  diagonal = Bind.join
  {-# INLINE diagonal #-}

instance Ord k => Trace (Map k) where
  diagonal = Bind.join
  {-# INLINE diagonal #-}

instance (Eq k, Hashable k) => Trace (HashMap k) where
  diagonal = Bind.join
  {-# INLINE diagonal #-}

instance Dim n => Trace (V n)
instance Trace V0
instance Trace V1
instance Trace V2
instance Trace V3
instance Trace V4
instance Trace Plucker
instance Trace Quaternion

instance Trace Complex where
  trace ((a :+ _) :+ (_ :+ b)) = a + b
  {-# INLINE trace #-}
  diagonal ((a :+ _) :+ (_ :+ b)) = a :+ b
  {-# INLINE diagonal #-}

instance (Trace f, Trace g) => Trace (Product f g) where
  trace (Pair xx yy) = trace (pfst <$> xx) + trace (psnd <$> yy) where
    pfst (Pair x _) = x
    psnd (Pair _ y) = y
  {-# INLINE trace #-}
  diagonal (Pair xx yy) = diagonal (pfst <$> xx) `Pair` diagonal (psnd <$> yy) where
    pfst (Pair x _) = x
    psnd (Pair _ y) = y
  {-# INLINE diagonal #-}

instance (Distributive g, Trace g, Trace f) => Trace (Compose g f) where
  trace = trace . fmap (fmap trace . distribute) . getCompose . fmap getCompose
  {-# INLINE trace #-}
  diagonal = Compose . fmap diagonal . diagonal . fmap distribute . getCompose . fmap getCompose
  {-# INLINE diagonal #-}

-- | Compute the <http://mathworld.wolfram.com/FrobeniusNorm.html Frobenius norm> of a matrix.
frobenius :: (Num a, Foldable f, Additive f, Additive g, Distributive g, Trace g) => f (g a) -> a
frobenius m = trace $ fmap (\ f' -> Foldable.foldl' (^+^) zero $ liftI2 (*^) f' m) (distribute m)

max_stars_repo_name: ekmett/linear

max_issues_repo_name: ekmett/linear

max_forks_repo_name: ekmett/linear


text:
module Histogram (
  histogram
) where

import qualified Data.Set as Set
import qualified Data.Vector.Generic as VG
import           Data.Vector.Unboxed ((!))
import qualified Statistics.Sample.Histogram as SSH
import           Statistics.Test.ApproxRand
import           Statistics.Types (Sample)

histogram ::
     Int
  -> TestResult
  -> Either String [(Double, Int)]
histogram prefBins (TestResult _ _ randomizedStats)
  | bins < 2  = Left "Cannot make a histogram - data cannot be divided in more than one bin."
  | otherwise = Right $ VG.toList $ VG.zip ticks numSamples
  where
    bins = numBins prefBins randomizedStats
    (lowerBounds, numSamples) =
      SSH.histogram prefBins randomizedStats

    -- Compute the middle values of the bins.
    bucketSize = (lowerBounds ! 1) - (lowerBounds ! 0)
    bucketHalf = bucketSize / 2
    ticks = VG.map ((+) bucketHalf) lowerBounds

numBins :: Int -> Sample -> Int
numBins pref =
  min pref . Set.size . VG.foldl' (flip Set.insert) Set.empty

max_stars_repo_name: danieldk/approx-rand-test

max_issues_repo_name: danieldk/approx-rand-test

max_forks_repo_name: danieldk/approx-rand-test


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ViewPatterns #-}
{- |
Module      : Language.Scheme.Types
Copyright   : Justin Ethier
Licence     : MIT (see LICENSE in the distribution)

Maintainer  : github.com/justinethier
Stability   : experimental
Portability : portable

This module contains top-level data type definitions, environments, error types, and associated functions.

-}

module Language.Scheme.Types
    (
    -- * Environments
      Env (..)
    , nullEnv
    -- * Error Handling
    , LispError (..)
    , ThrowsError
    , IOThrowsError
    , liftThrows
    , showCallHistory
    -- * Types and related functions
    , LispVal (
          Atom
        , List
        , DottedList
        , Vector
        , ByteVector
        , HashTable
        , Number
        , Float
        , Complex
        , Rational
        , String
        , Char
        , Bool
        , PrimitiveFunc
        , Func
             , params
             , vararg
             , body
             , closure
        , HFunc
             , hparams
             , hvararg
             , hbody
             , hclosure
        , IOFunc
        , CustFunc
        , EvalFunc
        , Pointer
             , pointerVar
             , pointerEnv
        , Opaque
        , Port
        , Continuation
             , contClosure
             , currentCont
             , nextCont
             , dynamicWind
             , contCallHist
        , Syntax
             , synClosure
             , synRenameClosure
             , synDefinedInMacro
             , synEllipsis
             , synIdentifiers
             , synRules
        , SyntaxExplicitRenaming
        , LispEnv
        , EOF
        , Nil)
    , nullLisp
    , toOpaque
    , fromOpaque
    , DeferredCode (..)
    , DynamicWinders (..)
    , makeNullContinuation
    , makeCPS
    , makeCPSWArgs
    , eqv
    , eqvList
    , eqVal
    , box
    , makeFunc
    , makeNormalFunc
    , makeVarargs
    , makeHFunc
    , makeNormalHFunc
    , makeHVarargs
    , validateFuncParams
    )
 where
import Control.Monad.Error
import Data.Complex
import Data.Array
import qualified Data.ByteString as BS
import Data.Dynamic
import qualified Data.Knob as DK
import qualified Data.List as DL
import Data.IORef
import qualified Data.Map
-- import Data.Maybe
import Data.Ratio
import System.IO
import Text.ParserCombinators.Parsec hiding (spaces)

-- Environment management

-- |A Scheme environment containing variable bindings of form @(namespaceName, variableName), variableValue@
data Env = Environment {
        parentEnv :: (Maybe Env),
        bindings :: (IORef (Data.Map.Map String (IORef LispVal))),
        pointers :: (IORef (Data.Map.Map String (IORef [LispVal])))
    }

instance Eq Env where
    (Environment _ xb xpts) == (Environment _ yb ypts) =
      (xb == yb) && (xpts == ypts)

-- |An empty environment
nullEnv :: IO Env
nullEnv = do
    nullBindings <- newIORef $ Data.Map.fromList []
    nullPointers <- newIORef $ Data.Map.fromList []
    return $ Environment Nothing nullBindings nullPointers

-- |Types of errors that may occur when evaluating Scheme code
data LispError = NumArgs (Maybe Integer) [LispVal] -- ^Invalid number of function arguments
  | TypeMismatch String LispVal -- ^Type error
  | Parser ParseError -- ^Parsing error
  | BadSpecialForm String LispVal -- ^Invalid special (built-in) form
  | UnboundVar String String -- ^ A referenced variable has not been declared
  | DivideByZero -- ^Divide by Zero error
  | NotImplemented String -- ^ Feature is not implemented
  | InternalError String {- ^An internal error within husk; in theory user (Scheme) code
                         should never allow one of these errors to be triggered. -}
  | Default String -- ^Default error
  | ErrorWithCallHist LispError [LispVal] -- ^Wraps an error to also include the current call history

-- |Create a textual description for a 'LispError'
showError :: LispError -> String
showError (NumArgs (Just expected) found) = "Expected " ++ show expected
                                  ++ " args but found "
                                  ++ (show $ length found)
                                  ++ " values: " ++ unwordsList found
showError (NumArgs Nothing found) = "Incorrect number of args, "
                                    ++ " found "
                                    ++ (show $ length found)
                                    ++ " values: " ++ unwordsList found
showError (TypeMismatch expected found) = "Invalid type: expected " ++ expected
                                  ++ ", found " ++ show found
showError (Parser parseErr) = "Parse error at " ++ ": " ++ show parseErr
showError (BadSpecialForm message form) = message ++ ": " ++ show form
showError (UnboundVar message varname) = message ++ ": " ++ varname
showError (DivideByZero) = "Division by zero"
showError (NotImplemented message) = "Not implemented: " ++ message
showError (InternalError message) = "An internal error occurred: " ++ message
showError (Default message) = "Error: " ++ message
showError (ErrorWithCallHist err stack) = showCallHistory (show err) stack

instance Show LispError where show = showError
instance Error LispError where
  noMsg = Default "An error has occurred"
  strMsg = Default

-- |Display call history for an error
showCallHistory :: String -> [LispVal] -> String
showCallHistory message hist = do
  let nums :: [Int]
      nums = [0..]
      ns = take (length hist) nums
  message ++ "\n\nCall History:\n" ++
    (unlines $ map (\(n, s) -> ('#' : show n) ++ ": " ++ show s)
                   (zip ns $ reverse hist))

-- |Container used by operations that could throw an error
type ThrowsError = Either LispError

-- |Container used to provide error handling in the IO monad
type IOThrowsError = ErrorT LispError IO

-- |Lift a ThrowsError into the IO monad
liftThrows :: ThrowsError a -> IOThrowsError a
liftThrows (Left err) = throwError err
liftThrows (Right val) = return val

-- |Scheme data types
data LispVal = Atom String
 -- ^Symbol
 | List [LispVal]
 -- ^List
 | DottedList [LispVal] LispVal
 -- ^Pair
 | Vector (Array Int LispVal)
 -- ^Vector
 | ByteVector BS.ByteString
 -- ^ByteVector from R7RS
 | HashTable (Data.Map.Map LispVal LispVal)
 {- ^Hash table.
 Technically this could be a derived data type instead of being built-in to the
 interpreter. And perhaps in the future it will be. But for now, a hash table
 is too important of a data type to not be included. -}
 --
 -- Map is technically the wrong structure to use for a hash table since it is based on a binary tree and hence operations tend to be O(log n) instead of O(1). However, according to <http://www.opensubscriber.com/message/haskell-cafe@haskell.org/10779624.html> Map has good performance characteristics compared to the alternatives. So it stays for the moment...
 --
 | Number Integer -- ^Integer number
 {- FUTURE: rename this to @Integer@ (or @WholeNumber@ or something else more meaningful)
 Integer -}
 | Float Double -- ^Double-precision floating point number
 {- FUTURE: rename this @Real@ instead of @Float@...
 Floating point -}
 | Complex (Complex Double)
 -- ^Complex number
 | Rational Rational
 -- ^Rational number
 | String String
 -- ^String
 | Char Char
 -- ^Character
 | Bool Bool
 -- ^Boolean
 | PrimitiveFunc ([LispVal] -> ThrowsError LispVal)
 -- ^Primitive function
 | Func {params :: [String],
         vararg :: (Maybe String),
         body :: [LispVal],
         closure :: Env
        }
 -- ^Function written in Scheme
 | HFunc {hparams :: [String],
          hvararg :: (Maybe String),
          hbody :: (Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal),
          hclosure :: Env
        }
 -- ^Function formed from a Haskell function
 | IOFunc ([LispVal] -> IOThrowsError LispVal)
 -- ^Primitive function within the IO monad
 | EvalFunc ([LispVal] -> IOThrowsError LispVal)
 {- ^Function within the IO monad with access to
 the current environment and continuation. -}
 | CustFunc ([LispVal] -> IOThrowsError LispVal)
 -- ^A custom function written by code outside of husk.
 --  Any code that uses the Haskell API should define custom
 --  functions using this data type.
 | Pointer { pointerVar :: String
            ,pointerEnv :: Env }
 -- ^Pointer to an environment variable.
 | Opaque Dynamic
 -- ^Opaque Haskell value.
 | Port Handle (Maybe DK.Knob)
 -- ^I/O port
 | Continuation {  contClosure :: Env                   -- Environment of the continuation
                 , currentCont :: (Maybe DeferredCode)  -- Code of current continuation
                 , nextCont :: (Maybe LispVal)          -- Code to resume after body of cont
                 , dynamicWind :: (Maybe [DynamicWinders]) -- Functions injected by (dynamic-wind)
                 , contCallHist :: [LispVal] -- Active call history
                }
 -- ^Continuation
 | Syntax { synClosure :: Maybe Env       -- ^ Code env in effect at definition time, if applicable
          , synRenameClosure :: Maybe Env -- ^ Renames (from macro hygiene) in effect at def time;
                                          --   only applicable if this macro defined inside another macro.
          , synDefinedInMacro :: Bool     -- ^ Set if macro is defined within another macro
          , synEllipsis :: String         -- ^ String to use as the ellipsis identifier
          , synIdentifiers :: [LispVal]   -- ^ Literal identifiers from syntax-rules
          , synRules :: [LispVal]         -- ^ Rules from syntax-rules
   } -- ^ Type to hold a syntax object that is created by a macro definition.
     --   Syntax objects are not used like regular types in that they are not
     --   passed around within variables. In other words, you cannot use set! to
     --   assign a variable to a syntax object. But they are used during function
     --   application. In any case, it is convenient to define the type here
     --   because syntax objects are stored in the same environments and
     --   manipulated by the same functions as regular variables.
 | SyntaxExplicitRenaming LispVal
   -- ^ Syntax for an explicit-renaming macro
 | LispEnv Env
   -- ^ Wrapper for a scheme environment
 | EOF
   -- ^ End of file indicator
 | Nil String
 -- ^Internal use only; do not use this type directly.

-- | Scheme /null/ value
nullLisp :: LispVal
nullLisp = List []

-- |Convert a Haskell value to an opaque Lisp value.
toOpaque :: Typeable a => a -> LispVal
toOpaque = Opaque . toDyn

-- |Convert an opaque Lisp value back into a Haskell value of the appropriate
--  type, or produce a TypeMismatch error.
fromOpaque :: forall a. Typeable a => LispVal -> ThrowsError a
-- fromOpaque (Opaque o) | isJust $ fromDynamic o = fromJust $ fromDynamic o
-- fromOpaque badArg = throwError $ TypeMismatch (show $ toOpaque (undefined :: a)) badArg

-- Old version that used ViewPatterns
fromOpaque (Opaque (fromDynamic -> Just v)) = return v
fromOpaque badArg = throwError $ TypeMismatch (show $ toOpaque (undefined :: a)) badArg

-- |Container to hold code that is passed to a continuation for deferred execution
data DeferredCode =
    SchemeBody [LispVal] | -- ^A block of Scheme code
    HaskellBody {
       contFunction :: (Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal)
     , contFunctionArgs :: (Maybe [LispVal]) -- Arguments to the higher-order function
    } -- ^A Haskell function

-- |Container to store information from a dynamic-wind
data DynamicWinders = DynamicWinders {
    before :: LispVal -- ^Function to execute when resuming continuation within extent of dynamic-wind
  , after :: LispVal -- ^Function to execute when leaving extent of dynamic-wind
}

showDWVal :: DynamicWinders -> String
showDWVal (DynamicWinders b a) = "(" ++ (show b) ++ " . " ++ (show a) ++ ")"

instance Show DynamicWinders where show = showDWVal

-- |Make an /empty/ continuation that does not contain any code
makeNullContinuation :: Env -> LispVal
makeNullContinuation env = Continuation env Nothing Nothing Nothing []

-- |Make a continuation that takes a higher-order function (written in Haskell)
makeCPS :: Env
        -- ^ Environment
        -> LispVal
        -- ^ Current continuation
        -> (Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal)
        -- ^ Haskell function
        -> LispVal
        -- ^ The Haskell function packaged as a LispVal
makeCPS env cont@(Continuation {contCallHist=hist}) cps = Continuation env (Just (HaskellBody cps Nothing)) (Just cont) (dynamicWind cont) hist
makeCPS env cont cps = Continuation env (Just (HaskellBody cps Nothing)) (Just cont) Nothing [] -- This overload just here for completeness; it should never be used

-- |Make a continuation that stores a higher-order function and arguments to that function
makeCPSWArgs :: Env
        -- ^ Environment
        -> LispVal
        -- ^ Current continuation
        -> (Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal)
        -- ^ Haskell function
        -> [LispVal]
        -- ^ Arguments to the function
        -> LispVal
        -- ^ The Haskell function packaged as a LispVal
makeCPSWArgs env cont@(Continuation {dynamicWind=dynWind,contCallHist=hist}) cps args =
    Continuation
        env
        (Just (HaskellBody cps (Just args)))
        (Just cont) dynWind hist
makeCPSWArgs env cont cps args =
    -- This overload just here for completeness; it should never be used
    Continuation
        env
        (Just (HaskellBody cps (Just args)))
        (Just cont) Nothing []

instance Ord LispVal where
  compare (Bool a) (Bool b) = compare a b
  compare (Number a) (Number b) = compare a b
  compare (Rational a) (Rational b) = compare a b
  compare (Float a) (Float b) = compare a b
  compare (String a) (String b) = compare a b
  compare (Char a) (Char b) = compare a b
  compare (Atom a) (Atom b) = compare a b
{- compare (DottedList xs x) (DottedList xs x) = compare a b
Vector
HashTable
List
Func
Others? -}
  compare a b = compare (show a) (show b) -- Hack (??): sort alphabetically when types differ or have no handlers

-- |Compare two 'LispVal' instances
eqv :: [LispVal]
    -- ^ A list containing two values to compare
    -> ThrowsError LispVal
    -- ^ Result wrapped as a Bool
eqv [(Bool arg1), (Bool arg2)] = return $ Bool $ arg1 == arg2
eqv [(Number arg1), (Number arg2)] = return $ Bool $ arg1 == arg2
eqv [(Complex arg1), (Complex arg2)] = return $ Bool $ arg1 == arg2
eqv [(Rational arg1), (Rational arg2)] = return $ Bool $ arg1 == arg2
eqv [(Float arg1), (Float arg2)] = return $ Bool $ arg1 == arg2
eqv [(String arg1), (String arg2)] = return $ Bool $ arg1 == arg2
eqv [(Char arg1), (Char arg2)] = return $ Bool $ arg1 == arg2
eqv [(Atom arg1), (Atom arg2)] = return $ Bool $ arg1 == arg2
eqv [(DottedList xs x), (DottedList ys y)] = eqv [List $ xs ++ [x], List $ ys ++ [y]]
eqv [(Vector arg1), (Vector arg2)] = eqv [List (elems arg1), List (elems arg2)]
eqv [(ByteVector a), (ByteVector b)] = return $ Bool $ a == b
eqv [(HashTable arg1), (HashTable arg2)] =
  eqv [List (map (\ (x, y) -> List [x, y]) $ Data.Map.toAscList arg1),
       List (map (\ (x, y) -> List [x, y]) $ Data.Map.toAscList arg2)]
--
-- This comparison function may be too simplistic. Basically we check to see if
-- functions have the same calling interface. If they do, then we compare the
-- function bodies for equality.
--
--FUTURE:
--
-- The real solution for this and many of the other comparison functions is to
-- assign memory locations to data. Then we can just compare memory locations
-- in cases such as this one. But that is a much larger change.
eqv [x@(Func _ _ xBody _), y@(Func _ _ yBody _)] = do
  if (show x) /= (show y)
     then return $ Bool False
     else eqvList eqv [List xBody, List yBody]
eqv [x@(HFunc{}), y@(HFunc{})] = do
  if (show x) /= (show y)
     then return $ Bool False
     else return $ Bool True
--
eqv [x@(PrimitiveFunc _), y@(PrimitiveFunc _)] = return $ Bool $ (show x) == (show y)
eqv [x@(IOFunc _), y@(IOFunc _)] = return $ Bool $ (show x) == (show y)
eqv [x@(CustFunc _), y@(CustFunc _)] = return $ Bool $ (show x) == (show y)
eqv [x@(EvalFunc _), y@(EvalFunc _)] = return $ Bool $ (show x) == (show y)
-- FUTURE: comparison of two continuations
eqv [l1@(List _), l2@(List _)] = eqvList eqv [l1, l2]
eqv [_, _] = return $ Bool False
eqv badArgList = throwError $ NumArgs (Just 2) badArgList

-- |Compare two lists of haskell values, using the given comparison function
eqvList :: ([LispVal] -> ThrowsError LispVal) -> [LispVal] -> ThrowsError LispVal
eqvList eqvFunc [(List arg1), (List arg2)] =
    return $ Bool $ (length arg1 == length arg2) &&
                    all eqvPair (zip arg1 arg2)
    where eqvPair (x1, x2) = case eqvFunc [x1, x2] of
                               Left _ -> False
                               Right (Bool val) -> val
                               _ -> False -- OK?
eqvList _ _ = throwError $ Default "Unexpected error in eqvList"

-- |A more convenient way to call /eqv/
eqVal :: LispVal -> LispVal -> Bool
eqVal a b = do
  let result = eqv [a, b]
  case result of
    Left _ -> False
    Right (Bool val) -> val
    _ -> False -- Is this OK?

instance Eq LispVal where
  x == y = eqVal x y

-- |Create a textual description of a 'LispVal'
showVal :: LispVal -> String
showVal (Nil _) = ""
showVal (EOF) = "#!EOF"
showVal (LispEnv _) = "<env>"
showVal (String contents) = "\"" ++ contents ++ "\""
showVal (Char chr) = [chr]
showVal (Atom name) = name
showVal (Number contents) = show contents
showVal (Complex contents) = show (realPart contents) ++ "+" ++ (show $ imagPart contents) ++ "i"
showVal (Rational contents) = (show (numerator contents)) ++ "/" ++ (show (denominator contents))
showVal (Float contents) = show contents
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (Vector contents) = "#(" ++ unwordsList (Data.Array.elems contents) ++ ")"
showVal (ByteVector contents) = "#u8(" ++ unwords (map show (BS.unpack contents)) ++ ")"
showVal (HashTable _) = "<hash-table>"
showVal (List contents) = "(" ++ unwordsList contents ++ ")"
showVal (DottedList h t) = "(" ++ unwordsList h ++ " . " ++ showVal t ++ ")"
showVal (PrimitiveFunc _) = "<primitive>"
showVal (Continuation {}) = "<continuation>"
showVal (Syntax {}) = "<syntax>"
showVal (SyntaxExplicitRenaming _) = "<er-macro-transformer syntax>"
showVal (Func {params = args, vararg = varargs, body = _, closure = _}) =
  "(lambda (" ++ unwords args ++
    (case varargs of
      Nothing -> ""
      Just arg -> " . " ++ arg) ++ ") ...)"
showVal (HFunc {hparams = args, hvararg = varargs, hbody = _, hclosure = _}) =
  "(lambda (" ++ unwords args ++
    (case varargs of
      Nothing -> ""
      Just arg -> " . " ++ arg) ++ ") ...)"
showVal (Port _ _) = "<IO port>"
showVal (IOFunc _) = "<IO primitive>"
showVal (CustFunc _) = "<custom primitive>"
showVal (EvalFunc _) = "<procedure>"
showVal (Pointer p _) = "<ptr " ++ p ++ ">"
showVal (Opaque d) = "<Haskell " ++ show (dynTypeRep d) ++ ">"

-- |A helper function to make pointer deref code more concise
box :: LispVal -> IOThrowsError [LispVal]
box a = return [a]

-- |Convert a list of Lisp objects into a space-separated string
unwordsList :: [LispVal] -> String
unwordsList = unwords . map showVal

-- |Allow conversion of lispval instances to strings
instance Show LispVal where show = showVal

-- Functions required by the interpreter --

-- |Create a scheme function
makeFunc :: -- forall (m :: * -> *).
            (Monad m) =>
            Maybe String -> Env -> [LispVal] -> [LispVal] -> m LispVal
makeFunc varargs env fparams fbody = return $ Func (map showVal fparams) varargs fbody env

-- |Create a normal scheme function
makeNormalFunc :: (Monad m) => Env
               -> [LispVal]
               -> [LispVal]
               -> m LispVal
makeNormalFunc = makeFunc Nothing

-- |Create a scheme function that can receive any number of arguments
makeVarargs :: (Monad m) => LispVal -> Env
                        -> [LispVal]
                        -> [LispVal]
                        -> m LispVal
makeVarargs = makeFunc . Just . showVal

-- Functions required by a compiled program --

-- |Create a haskell function
makeHFunc ::
            (Monad m) =>
            Maybe String
         -> Env
         -> [String]
         -> (Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal)
--         -> String
         -> m LispVal
makeHFunc varargs env fparams fbody = return $ HFunc fparams varargs fbody env --(map showVal fparams) varargs fbody env
-- |Create a normal haskell function
makeNormalHFunc :: (Monad m) =>
                  Env
               -> [String]
               -> (Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal)
               -> m LispVal
makeNormalHFunc = makeHFunc Nothing

-- |Create a haskell function that can receive any number of arguments
makeHVarargs :: (Monad m) => LispVal
                        -> Env
                        -> [String]
                        -> (Env -> LispVal -> LispVal -> Maybe [LispVal] -> IOThrowsError LispVal)
                        -> m LispVal
makeHVarargs = makeHFunc . Just . showVal

-- |Validate formal function parameters.
validateFuncParams :: [LispVal] -> Maybe Integer -> IOThrowsError Bool
--validateFuncParams [Atom _] _ = return True
validateFuncParams ps (Just n) = do
  if length ps /= fromInteger n
     then throwError $ NumArgs (Just n) ps
     else validateFuncParams ps Nothing
validateFuncParams ps Nothing = do
  let syms = filter filterArgs ps
  if (length syms) /= (length ps)
     then throwError $ Default $
             "Invalid lambda parameter(s): " ++ show (List ps)
     else do
         let strs = DL.sort $ map (\ (Atom a) -> a) ps
         case dupe strs of
            Just d -> throwError $ Default $
                         "Duplicate lambda parameter " ++ d
            _ -> return True
 where
  filterArgs (Atom _) = True
  filterArgs _ = False

  dupe (a : b : rest)
    | a == b = Just a
    | otherwise = dupe (b : rest)
  dupe _ = Nothing

max_stars_repo_name: farrellm/husk-scheme

max_issues_repo_name: haskell-lisp/husk-scheme

max_forks_repo_name: haskell-lisp/husk-scheme