text:
{-# LANGUAGE CPP, DerivingVia, GADTs, InstanceSigs, KindSignatures    #-}
{-# LANGUAGE PatternSynonyms, RankNTypes, RoleAnnotations             #-}
{-# LANGUAGE ScopedTypeVariables, StandaloneDeriving, TemplateHaskell #-}
{-# LANGUAGE TypeApplications, TypeFamilies, TypeOperators            #-}
{-# LANGUAGE UndecidableSuperClasses                                  #-}
{-# OPTIONS_GHC -Wno-orphans #-}
module Control.Subcategory.Functor
  ( Constrained(..), Dom(), CFunctor (..),
    (<$:>),
    defaultCmapConst,
    WrapFunctor (..),
    WrapMono (WrapMono, unwrapMono),
    coerceToMono, withMonoCoercible,
  )
where
import qualified Control.Applicative                  as App
import           Control.Arrow                        (Arrow, ArrowMonad)
import           Control.Exception                    (Handler)
import qualified Control.Monad.ST.Lazy                as LST
import qualified Control.Monad.ST.Strict              as SST
import           Control.Subcategory.Wrapper.Internal
import           Data.Coerce
import           Data.Complex                         (Complex)
import qualified Data.Functor.Compose                 as SOP
import           Data.Functor.Const                   (Const)
import           Data.Functor.Identity                (Identity)
import qualified Data.Functor.Product                 as SOP
import qualified Data.Functor.Sum                     as SOP
import           Data.Hashable                        (Hashable)
import qualified Data.HashMap.Strict                  as HM
import qualified Data.HashSet                         as HS
import qualified Data.IntMap                          as IM
import           Data.Kind                            (Constraint, Type)
import           Data.List.NonEmpty                   (NonEmpty)
import qualified Data.Map                             as Map
import qualified Data.Monoid                          as Mon
import           Data.MonoTraversable                 (Element,
                                                       MonoFunctor (..))
#if MIN_VERSION_mono_traversable(1,0,14)
import Data.MonoTraversable (WrappedMono)
#endif

import qualified Data.IntSet                     as IS
import           Data.Ord                        (Down (..))
import qualified Data.Primitive.Array            as A
import qualified Data.Primitive.PrimArray        as PA
import qualified Data.Primitive.SmallArray       as SA
import           Data.Proxy                      (Proxy)
import qualified Data.Semigroup                  as Sem
import qualified Data.Sequence                   as Seq
import qualified Data.Set                        as Set
import qualified Data.Tree                       as Tree
import qualified Data.Vector                     as V
import qualified Data.Vector.Primitive           as P
import qualified Data.Vector.Storable            as S
import qualified Data.Vector.Unboxed             as U
import           Foreign.Ptr                     (Ptr)
import           GHC.Conc                        (STM)
import           GHC.Generics                    ((:*:) (..), (:+:) (..),
                                                  (:.:) (..), K1, M1, Par1,
                                                  Rec1, U1, URec, V1)
import qualified System.Console.GetOpt           as GetOpt
import           Text.ParserCombinators.ReadP    (ReadP)
import           Text.ParserCombinators.ReadPrec (ReadPrec)

infixl 4 <$:

class Constrained (f :: Type -> Type) where
  type Dom f (a :: Type) :: Constraint
  type Dom f a = ()

class Constrained f => CFunctor f where
  cmap :: (Dom f a, Dom f b) => (a -> b) -> f a -> f b
  default cmap :: Functor f => (a -> b) -> f a -> f b
  cmap = fmap
  {-# INLINE cmap #-}
  (<$:) :: (Dom f a, Dom f b) => a -> f b -> f a
  (<$:) = cmap . const
  {-# INLINE (<$:) #-}

defaultCmapConst :: (CFunctor f, Dom f a, Dom f b) => a -> f b -> f a
defaultCmapConst = cmap . const
{-# INLINE defaultCmapConst #-}

instance Constrained (WrapFunctor f) where
  type Dom (WrapFunctor f) a = ()

instance Functor f => CFunctor (WrapFunctor f) where
  cmap :: (a -> b) -> WrapFunctor f a -> WrapFunctor f b
  cmap = fmap
  {-# INLINE cmap #-}
  (<$:) :: a -> WrapFunctor f b -> WrapFunctor f a
  (<$:) = (<$)
  {-# INLINE (<$:) #-}

instance Constrained []
instance CFunctor []
instance Constrained Maybe
instance CFunctor Maybe
instance Constrained IO
instance CFunctor IO
instance Constrained Par1
instance CFunctor Par1
instance Constrained NonEmpty
instance CFunctor NonEmpty
instance Constrained ReadP
instance CFunctor ReadP
instance Constrained ReadPrec
instance CFunctor ReadPrec

instance Constrained Down
instance CFunctor Down
instance Constrained Mon.Product
instance CFunctor Mon.Product

instance Constrained Mon.Sum
instance CFunctor Mon.Sum
instance Constrained Mon.Dual
instance CFunctor Mon.Dual

instance Constrained Mon.Last
instance CFunctor Mon.Last
instance Constrained Mon.First
instance CFunctor Mon.First

instance Constrained STM
instance CFunctor STM
instance Constrained Handler
instance CFunctor Handler

instance Constrained Identity
instance CFunctor Identity
instance Constrained App.ZipList
instance CFunctor App.ZipList
instance Constrained GetOpt.ArgDescr
instance CFunctor GetOpt.ArgDescr
instance Constrained GetOpt.OptDescr
instance CFunctor GetOpt.OptDescr
instance Constrained GetOpt.ArgOrder
instance CFunctor GetOpt.ArgOrder
instance Constrained Sem.Option
instance CFunctor Sem.Option

instance Constrained Sem.Last
instance CFunctor Sem.Last
instance Constrained Sem.First
instance CFunctor Sem.First

instance Constrained Sem.Max
instance CFunctor Sem.Max
instance Constrained Sem.Min
instance CFunctor Sem.Min

instance Constrained Complex
instance CFunctor Complex
instance Constrained (Either a)
instance CFunctor (Either a)

instance Constrained V1
instance CFunctor V1
instance Constrained U1
instance CFunctor U1

instance Constrained ((,) a)
instance CFunctor ((,) a)
instance Constrained (SST.ST s)
instance CFunctor (SST.ST s)

instance Constrained (LST.ST s)
instance CFunctor (LST.ST s)
instance Constrained Proxy
instance CFunctor Proxy

instance Constrained (ArrowMonad a)
instance Arrow a => CFunctor (ArrowMonad a)
instance Constrained (App.WrappedMonad m)
instance Monad m => CFunctor (App.WrappedMonad m)

instance Constrained (Sem.Arg a)
instance CFunctor (Sem.Arg a)
instance Constrained (Rec1 f)
instance Functor f => CFunctor (Rec1 f)

instance Constrained (URec Char)
instance CFunctor (URec Char)
instance Constrained (URec Double)
instance CFunctor (URec Double)

instance Constrained (URec Float)
instance CFunctor (URec Float)
instance Constrained (URec Int)
instance CFunctor (URec Int)

instance Constrained (URec Word)
instance CFunctor (URec Word)
instance Constrained (URec (Ptr ()))
instance CFunctor (URec (Ptr ()))

instance Constrained f => Constrained (Mon.Ap f) where
  type Dom (Mon.Ap f) a = Dom f a

deriving newtype instance CFunctor f => CFunctor (Mon.Ap f)

instance Constrained (Mon.Alt f) where
  type Dom (Mon.Alt f) a = Dom f a
deriving newtype instance CFunctor f => CFunctor (Mon.Alt f)

instance Constrained (Const m)
instance CFunctor (Const m)
instance Constrained (App.WrappedArrow a b)
instance Arrow a => CFunctor (App.WrappedArrow a b)

instance Constrained ((->) r)
instance CFunctor ((->) r)
instance Constrained (K1 i c)
instance CFunctor (K1 i c)

instance Constrained (f :+: g) where
  type Dom (f :+: g) a = (Dom f a, Dom g a)
instance (CFunctor f, CFunctor g) => CFunctor (f :+: g) where
  cmap f (L1 xs) = L1 $ cmap f xs
  cmap f (R1 xs) = R1 $ cmap f xs
  {-# INLINE [1] cmap #-}
instance Constrained (f :*: g) where
  type Dom (f :*: g) a = (Dom f a, Dom g a)
instance (CFunctor f, CFunctor g) => CFunctor (f :*: g) where
  cmap f (l :*: r) = cmap f l :*: cmap f r
  {-# INLINE cmap #-}

instance Constrained (f :.: (g :: Type -> Type)) where
  type Dom (f :.: g) a = (Dom f (g a), Dom g a)
instance (CFunctor f, CFunctor g) => CFunctor (f :.: g) where
  cmap f gfa = Comp1 $ cmap (cmap f) $ unComp1 gfa
  {-# INLINE cmap #-}
instance (Constrained f, Constrained g) => Constrained (SOP.Sum f g) where
  type Dom (SOP.Sum f g) a = (Dom f a, Dom g a)

instance (CFunctor f, CFunctor g) => CFunctor (SOP.Sum f g) where
  cmap f (SOP.InL a) = SOP.InL $ cmap f a
  cmap f (SOP.InR b) = SOP.InR $ cmap f b
  {-# INLINE cmap #-}

  (<$:) = defaultCmapConst
  {-# INLINE (<$:) #-}

instance (Constrained f, Constrained g) => Constrained (SOP.Product f g) where
  type Dom (SOP.Product f g) a = (Dom f a, Dom g a)

instance (CFunctor f, CFunctor g) => CFunctor (SOP.Product f g) where
  cmap f (SOP.Pair a b) = SOP.Pair (cmap f a) (cmap f b)
  {-# INLINE cmap #-}

  (<$:) = defaultCmapConst
  {-# INLINE (<$:) #-}

instance (Constrained (f ::Type -> Type), Constrained (g :: Type -> Type))
  => Constrained (SOP.Compose f g) where
  type Dom (SOP.Compose f g) a = (Dom g a, Dom f (g a))

instance (CFunctor f, CFunctor g) => CFunctor (SOP.Compose f g) where
  cmap f (SOP.Compose a) = SOP.Compose $ cmap (cmap f) a
  (<$:) = defaultCmapConst

  {-# INLINE (<$:) #-}

instance Constrained (M1 i c f)
instance Functor f => CFunctor (M1 i c f)

instance Constrained Seq.Seq
instance CFunctor Seq.Seq

#if MIN_VERSION_mono_traversable(1,0,14)
instance Constrained (WrappedMono mono) where
  type Dom (WrappedMono mono) a = a ~ Element mono

instance MonoFunctor IS.IntSet where
  omap = IS.map

instance MonoFunctor mono => CFunctor (WrappedMono mono) where
  cmap = omap
  (<$:) = omap . const
#endif

instance Constrained (WrapMono mono) where
  type Dom (WrapMono mono) b = b ~ Element mono

instance {-# OVERLAPPABLE #-} MonoFunctor a
      => CFunctor (WrapMono a) where
  cmap = coerce @((Element a -> Element a) -> a -> a) omap
  {-# INLINE [1] cmap #-}

  (<$:) = defaultCmapConst
  {-# INLINE [1] (<$:) #-}

instance Constrained IM.IntMap
instance CFunctor IM.IntMap

instance Constrained (Map.Map k)
instance Ord k => CFunctor (Map.Map k)

instance Constrained Set.Set where
  type Dom Set.Set a = Ord a

instance CFunctor Set.Set where
  cmap = Set.map
  {-# INLINE [1] cmap #-}
  (<$:) = flip $ \s ->
    if Set.null s
    then const Set.empty else Set.singleton
  {-# INLINE [1] (<$:) #-}

instance Constrained HS.HashSet where
  type Dom HS.HashSet a = (Hashable a, Eq a)

instance CFunctor HS.HashSet where
  cmap :: (Hashable b, Eq b) => (a -> b) -> HS.HashSet a -> HS.HashSet b
  cmap = HS.map
  {-# INLINE [1] cmap #-}
  (<$:) = flip $ \s -> if HS.null s
    then const HS.empty else HS.singleton
  {-# INLINE (<$:) #-}

instance Constrained (HM.HashMap k)
instance CFunctor (HM.HashMap k)
instance Constrained Tree.Tree
instance CFunctor Tree.Tree

infixl 4 <$:>
(<$:>) :: (CFunctor f, Dom f a, Dom f b) => (a -> b) -> f a -> f b
(<$:>) = cmap
{-# INLINE [1] (<$:>) #-}

instance Constrained V.Vector where
  type Dom V.Vector a = ()

instance CFunctor V.Vector where
  cmap = V.map
  {-# INLINE [1] cmap #-}

instance Constrained U.Vector where
  type Dom U.Vector a = U.Unbox a
instance CFunctor U.Vector where
  cmap = U.map
  {-# INLINE [1] cmap #-}
instance Constrained S.Vector where
  type Dom S.Vector a = S.Storable a
instance CFunctor S.Vector where
  cmap = S.map
  {-# INLINE [1] cmap #-}

instance Constrained P.Vector where
  type Dom P.Vector a = P.Prim a
instance CFunctor P.Vector where
  cmap = P.map
  {-# INLINE [1] cmap #-}

instance Constrained PA.PrimArray where
  type Dom PA.PrimArray a = P.Prim a

instance CFunctor PA.PrimArray where
  cmap = PA.mapPrimArray
  {-# INLINE [1] cmap #-}

deriving via WrapFunctor SA.SmallArray
  instance Constrained SA.SmallArray
deriving via WrapFunctor SA.SmallArray
  instance CFunctor SA.SmallArray

deriving via WrapFunctor A.Array
  instance Constrained A.Array
deriving via WrapFunctor A.Array
  instance CFunctor A.Array

max_stars_repo_name: konn/subcategories

max_issues_repo_name: konn/subcategories

max_forks_repo_name: konn/subcategories


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Statistics.Viterbi.Discrete.Lazy
-- Copyright   :  (c) A. V. H. McPhail 2011
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- Viterbi Algorithm
--
-----------------------------------------------------------------------------

module Numeric.Statistics.Viterbi.Discrete.Lazy (
                          ) where

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

import qualified Data.Array.IArray as I

import Numeric.LinearAlgebra

import Numeric.GSL.Statistics

import Numeric.Statistics

import System.Random

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

newtype Viterbi a = Viterbi {
      _memory :: Int
    , _states :: Int

train :: Int       -- ^ sequence length

ext: hs

max_stars_repo_name: amcphail/hstatistics

max_issues_repo_name: amcphail/hstatistics

max_forks_repo_name: amcphail/hstatistics


text:
{-# LANGUAGE BangPatterns #-}
module Game.Physics.Integration ( integrate
                                ) where

import Prelewd

import Control.Monad.ST
import Data.Packed.Development ( at' )
import Data.Packed.ST
import Data.Tuple
import Foreign.Storable ( Storable )
import qualified Data.Vector as D -- oh god this'll get confusing.
import Game.Physics.Types
import Numeric.GSL
import Numeric.LinearAlgebra hiding ( i )

-- | Converts a single physics object into its corresponding matrix
--   representation. This works for physics objects in any dimension.
--
--   > PObject (s_x s_y s_z)
--   >         (v_x v_y v_z)
--   >
--   > [ s_x v_x ]T
--   > | s_y v_y |
--   > [ s_z v_z ]
asMatrix :: PhysicsObject -> Matrix Double
asMatrix p = fromRows [ s, v ]
    where
        s = pos p
        v = vel p

-- | Converts a list of physics objects into their corresponding combined
--   matrix by stacking their 'asMatrix' representations vertically.
bigMat :: [PhysicsObject] -> Matrix Double
bigMat = fromColumns . concatMap (toColumns . asMatrix)

-- | This function might be in Haskell, but it's just as dangerous as C.
memcpy :: Storable t
       => STVector s t -- ^ dest
       -> Int -- ^ destination offset
       -> Vector t -- ^ src
       -> Int -- ^ source offset
       -> Int -- ^ number of elements to copy
       -> ST s ()
memcpy dst dstoff src srcoff n = for 0 (<n) (+1) $
                                    \idx -> unsafeWriteVector dst (dstoff+idx) (src `at'` (srcoff+idx))
{-# INLINE memcpy #-}

-- | An emulation of a C for loop.
--
--   for (int i = first clause; second clause; third clause)
--       { loop body }
for :: Monad m
    => a -- First clause.
    -> (a -> Bool) -- Second clause.
    -> (a -> a) -- Third clause.
    -> (a -> m b) -- Loop body. Return value is ignored for convenience.
    -> m ()
for i0 condition inc body = go i0
    where
        go !i | condition i = body i >> go (inc i)
              | otherwise   = return ()
        {-# INLINE go #-}
{-# INLINE for #-}

integrate :: D.Vector PhysicsObject -> Double -> D.Vector PhysicsObject
integrate ps _ | D.length ps == 0 = ps
integrate ps t = D.map (\((PObject _ _ im f), [s', v']) -> PObject s' v' im f)
               . D.zip ps
               . D.map toRows
               . D.fromList -- vector fusion can't handle concatenation.
               . concat
               . toBlocksEvery 2 dims
               . reshape width
               . (! 1) . toRows  -- second row.
               . odeSolveV RKf45  -- Is there a better choice of integrator?
                           t      -- Not really sure about this. odeSolve uses t/100.
                           1.0e-5 -- arbitrary. Lower it if simulation blows
                                  -- up due to inaccuracy. Highten it for
                                  -- inceased performance.
                           1.0e-5 -- same as above.
                           integrator
                           Nothing -- TODO: Can we supply a jacobian one day?
                                   --       Is it recommended, or totally unnecessary?
                                   --       I have no idea. What the fuck is a jacobian?
                           initialConditions
                           $ fromList [0.0, t]
    where
        initialConditions :: Vector Double
        initialConditions = flatten . bigMat $ D.toList ps

        -- integrates by shifting all values over one to the left, and querying
        -- all the force functions for new accelerations. This is written
        -- imperatively since it rapes the GC if done functionally, and takes
        -- 2x as long as a result.
        --
        -- msrc matrix layout (after a (reshape width))
        --
        -- [ s_x_0 s_y_0 s_z_0   s_x_1 s_y_1 s_z_1  ... ]
        -- [ v_x_0 v_y_0 v_z_0   v_x_1 v_y_1 v_z_1  ... ]
        --
        integrator :: Double -> Vector Double -> Vector Double
        integrator dt vsrc = runSTVector $ do
                vdst <- newUndefinedVector (dim vsrc)
                memcpy vdst width vsrc 0 width

                let xs = (subVector 0     width vsrc)
                    vs = (subVector width width vsrc)
                    objectAt i = let j = i*dims
                                     p = ps D.! i
                                  in PObject (subVector j dims xs)
                                             (subVector j dims vs)
                                             (invMass p)
                                             (netForce p)

                for (0, 0) ((< len) . fst) ((+1) *** (+dims)) $ \(i, j) ->
                    let p  = ps D.! i
                        im = invMass p
                        f  = netForce p objectAt dt
                     in for 0 (<dims) (+1) $ \k ->
                            unsafeWriteVector vdst (j+k) (im * (f `at'` k))

                return vdst

        -- dimensions of each point/vector.
        dims = dim . pos $ ps D.! 0

        -- the number of particles we're integrating.
        len = D.length ps

        -- width of every row in the matrix.
        width = dims*len

max_stars_repo_name: bfops/Chess

max_issues_repo_name: bfops/Chess

max_forks_repo_name: bfops/Chess


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

---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------
-- |
-- | Module : Benchmark approx median
-- | Creator: Xiao Ling
-- | Created: 12/14/2015
-- | see    : http://www.serpentine.com/criterion/tutorial.html
-- | TODO   : move this to benchmark directory
-- |
---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------

module BapproxMedian where

import Data.List
import Data.Random
import Data.Conduit
import qualified Data.Conduit.List as Cl
import Criterion.Main

import Core
import Statistics
import ApproxMedian

main :: IO ()
main = return ()


text:
{- |
This module provides  normalized multi-dimensional versions of the transforms in @fftw@.

The forwards transforms in this module are identical to those in "Numeric.FFT.Vector.Unnormalized".
The backwards transforms are normalized to be their inverse operations (approximately, due to floating point precision).

For more information on the underlying transforms, see
<http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html>.

@since 0.2
-}

module Numeric.FFT.Vector.Invertible.Multi
  (
        -- * Creating and executing 'Plan's
        run,
        plan,
        execute,
        -- * Complex-to-complex transforms
        U.dft,
        idft,
        -- * Real-to-complex transforms
        U.dftR2C,
        dftC2R,
  ) where

import Numeric.FFT.Vector.Base
import qualified Numeric.FFT.Vector.Unnormalized.Multi as U
import Data.Complex
import qualified Data.Vector.Storable as VS

-- | A backward discrete Fourier transform which is the inverse of 'U.dft'.  The output and input sizes are the same (@n@).
idft :: TransformND (Complex Double) (Complex Double)
idft = U.idft {normalizationND = \ns -> constMultOutput $ 1 / toEnum (VS.product ns)}

-- | A normalized backward discrete Fourier transform which is the left inverse of
-- 'U.dftR2C'.  (Specifically, @run dftC2R . run dftR2C == id@.)
--
-- This 'Transform' behaves differently than the others:
--
--  - Calling @planND dftC2R dims@ where @dims = [n0, ..., nk]@ creates a 'Plan' whose /output/ size is @dims@, and whose
--    /input/ size is @[n0, ..., nk \`div\` 2 + 1]@.
--
--  - If @length v == n0 * ... * nk@, then @length (run dftC2R v) == n0 * ... * 2*(nk-1)@.
--
dftC2R :: TransformND (Complex Double) Double
dftC2R = U.dftC2R {normalizationND = \ns -> constMultOutput $ 1 / toEnum (VS.product ns)}

max_stars_repo_name: TravisWhitaker/vector-fftw

max_issues_repo_name: TravisWhitaker/vector-fftw

max_forks_repo_name: TravisWhitaker/vector-fftw


text:
import Data.Complex
import Data.Ratio

-- types

a1 :: Int
a1 = 4

a2 :: Complex Double
a2 = 2 :+ 3

a3 :: Ratio Int
a3 = 2%3

main = do
       print a1
       print a2
       print a3

-- string literal

c1 :: Char
c1 = '\x0053'

c2 :: [Char]
c2 = "\o0120"

main2 = do
        putChar c1 >> putChar '\n'
        putChar (c2!!0) >> putChar '\n'

-- lists

l0 :: [Int]
l0 = [1, 3..10]

l1 :: [Int]
l1 = [1..]

l2 :: [String]
l2 = ["aa", "bb", "cc"]

main3 = do
       print l0
       print $ take 20 l1
       print l2

-- tuples

t0 :: (Int, Int)
t0 = (69,96)

t1 :: (String, Int, Double)
t1  = ("something", 8, 3.123)

t2 :: ([Int],([String], (Int, Float)))
t2 = ([1..4], (["1..4","sdf"], (1234, 1234.1234)))

main4 = do
        print t0
        print t1
        print t2

--functions

f0 :: String -> Int
f0 = length

f1 :: String -> (String, Int)
f1 x = (x, length x)

f2 :: [String] -> [(String, Int)]
f2 = map f1

main5 = do
        print $ f0 "1234567890"
        print $ f1 "0987654321"
        print $ f2 ["HEYYY", "Maaaan"]


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeSynonymInstances #-}

#if __GLASGOW_HASKELL__ >= 701
{-# LANGUAGE DefaultSignatures #-}
{-# LANGUAGE Trustworthy #-}
#endif

#if __GLASGOW_HASKELL__ >= 708
{-# LANGUAGE EmptyCase #-}
#endif

module Generics.Deriving.Show (
  -- * Generic show class
    GShow(..)

  -- * Default definition
  , gshowsPrecdefault

  -- * Internal show class
  , GShow'(..)

  ) where

import           Control.Applicative (Const, ZipList)

import           Data.Char (GeneralCategory)
import           Data.Int
import           Data.Monoid (All, Any, Dual, Product, Sum)
import qualified Data.Monoid as Monoid (First, Last)
import           Data.Version (Version)
import           Data.Word

import           Foreign.C.Types
import           Foreign.ForeignPtr (ForeignPtr)
import           Foreign.Ptr

import           Generics.Linear

import           GHC.Exts hiding (Any)

import           System.Exit (ExitCode)
import           System.IO (BufferMode, Handle, HandlePosn, IOMode, SeekMode)
import           System.IO.Error (IOErrorType)
import           System.Posix.Types

#if MIN_VERSION_base(4,4,0)
import           Data.Complex (Complex)
#endif

#if MIN_VERSION_base(4,7,0)
import           Data.Proxy (Proxy)
#endif

#if MIN_VERSION_base(4,8,0)
import           Data.Functor.Identity (Identity)
import           Data.Monoid (Alt)
import           Data.Void (Void)
import           Numeric.Natural (Natural)
#endif

#if MIN_VERSION_base(4,9,0)
import           Data.List.NonEmpty (NonEmpty)
import qualified Data.Semigroup as Semigroup (First, Last)
import           Data.Semigroup (Arg, Max, Min, WrappedMonoid)
#endif

--------------------------------------------------------------------------------
-- Generic show
--------------------------------------------------------------------------------

intersperse :: a -> [a] -> [a]
intersperse _ []    = []
intersperse _ [h]   = [h]
intersperse x (h:t) = h : x : (intersperse x t)

appPrec :: Int
appPrec = 2

data Type = Rec | Tup | Pref | Inf String

class GShow' f where
  gshowsPrec' :: Type -> Int -> f a -> ShowS
  isNullary   :: f a -> Bool
  isNullary = error "generic show (isNullary): unnecessary case"

instance GShow' V1 where
  gshowsPrec' _ _ x = case x of
#if __GLASGOW_HASKELL__ >= 708
                        {}
#else
                        !_ -> error "Void gshowsPrec"
#endif

instance GShow' U1 where
  gshowsPrec' _ _ U1 = id
  isNullary _ = True

instance (GShow c) => GShow' (K1 i c) where
  gshowsPrec' _ n (K1 a) = gshowsPrec n a
  isNullary _ = False

-- No instances for P or Rec because gshow is only applicable to types of kind *

instance (GShow' a, Constructor c) => GShow' (M1 C c a) where
  gshowsPrec' _ n c@(M1 x) =
    case fixity of
      Prefix    -> showParen (n > appPrec && not (isNullary x))
                    ( showString (conName c)
                    . if (isNullary x) then id else showChar ' '
                    . showBraces t (gshowsPrec' t appPrec x))
      Infix _ m -> showParen (n > m) (showBraces t (gshowsPrec' t m x))
      where fixity = conFixity c
            t = if (conIsRecord c) then Rec else
                  case (conIsTuple c) of
                    True -> Tup
                    False -> case fixity of
                                Prefix    -> Pref
                                Infix _ _ -> Inf (show (conName c))
            showBraces :: Type -> ShowS -> ShowS
            showBraces Rec     p = showChar '{' . p . showChar '}'
            showBraces Tup     p = showChar '(' . p . showChar ')'
            showBraces Pref    p = p
            showBraces (Inf _) p = p

            conIsTuple :: C1 c f p -> Bool
            conIsTuple y = tupleName (conName y) where
              tupleName ('(':',':_) = True
              tupleName _           = False

instance (Selector s, GShow' a) => GShow' (M1 S s a) where
  gshowsPrec' t n s@(M1 x) | selName s == "" = --showParen (n > appPrec)
                                                 (gshowsPrec' t n x)
                           | otherwise       =   showString (selName s)
                                               . showString " = "
                                               . gshowsPrec' t 0 x
  isNullary (M1 x) = isNullary x

instance GShow' a => GShow' (M1 D d a) where
  gshowsPrec' t n (M1 x) = gshowsPrec' t n x

instance GShow' a => GShow' (MP1 m a) where
  gshowsPrec' t n (MP1 x) = gshowsPrec' t n x

instance (GShow' a, GShow' b) => GShow' (a :+: b) where
  gshowsPrec' t n (L1 x) = gshowsPrec' t n x
  gshowsPrec' t n (R1 x) = gshowsPrec' t n x

instance (GShow' a, GShow' b) => GShow' (a :*: b) where
  gshowsPrec' t@Rec     n (a :*: b) =
    gshowsPrec' t n     a . showString ", " . gshowsPrec' t n     b
  gshowsPrec' t@(Inf s) n (a :*: b) =
    gshowsPrec' t n     a . showString s    . gshowsPrec' t n     b
  gshowsPrec' t@Tup     n (a :*: b) =
    gshowsPrec' t n     a . showChar ','    . gshowsPrec' t n     b
  gshowsPrec' t@Pref    n (a :*: b) =
    gshowsPrec' t (n+1) a . showChar ' '    . gshowsPrec' t (n+1) b

  -- If we have a product then it is not a nullary constructor
  isNullary _ = False

-- Unboxed types
instance GShow' UChar where
  gshowsPrec' _ _ (UChar c)   = showsPrec 0 (C# c) . showChar '#'
instance GShow' UDouble where
  gshowsPrec' _ _ (UDouble d) = showsPrec 0 (D# d) . showString "##"
instance GShow' UFloat where
  gshowsPrec' _ _ (UFloat f)  = showsPrec 0 (F# f) . showChar '#'
instance GShow' UInt where
  gshowsPrec' _ _ (UInt i)    = showsPrec 0 (I# i) . showChar '#'
instance GShow' UWord where
  gshowsPrec' _ _ (UWord w)   = showsPrec 0 (W# w) . showString "##"

class GShow a where
  gshowsPrec :: Int -> a -> ShowS
#if __GLASGOW_HASKELL__ >= 701
  default gshowsPrec :: (Generic a, GShow' (Rep a))
                     => Int -> a -> ShowS
  gshowsPrec = gshowsPrecdefault
#endif

  gshows :: a -> ShowS
  gshows = gshowsPrec 0

  gshow :: a -> String
  gshow x = gshows x ""

  gshowList :: [a] -> ShowS
  gshowList l =   showChar '['
                . foldr (.) id
                   (intersperse (showChar ',') (map (gshowsPrec 0) l))
                . showChar ']'

gshowsPrecdefault :: (Generic a, GShow' (Rep a))
                  => Int -> a -> ShowS
gshowsPrecdefault n = gshowsPrec' Pref n . from

-- Base types instances
-- Base types instances
instance GShow () where
  gshowsPrec = gshowsPrecdefault

instance (GShow a, GShow b) => GShow (a, b) where
  gshowsPrec = gshowsPrecdefault

instance (GShow a, GShow b, GShow c) => GShow (a, b, c) where
  gshowsPrec = gshowsPrecdefault

instance (GShow a, GShow b, GShow c, GShow d) => GShow (a, b, c, d) where
  gshowsPrec = gshowsPrecdefault

instance (GShow a, GShow b, GShow c, GShow d, GShow e) => GShow (a, b, c, d, e) where
  gshowsPrec = gshowsPrecdefault

instance (GShow a, GShow b, GShow c, GShow d, GShow e, GShow f)
    => GShow (a, b, c, d, e, f) where
  gshowsPrec = gshowsPrecdefault

instance (GShow a, GShow b, GShow c, GShow d, GShow e, GShow f, GShow g)
    => GShow (a, b, c, d, e, f, g) where
  gshowsPrec = gshowsPrecdefault

instance GShow a => GShow [a] where
  gshowsPrec _ = gshowList

instance (GShow (f p), GShow (g p)) => GShow ((f :+: g) p) where
  gshowsPrec = gshowsPrecdefault

instance (GShow (f p), GShow (g p)) => GShow ((f :*: g) p) where
  gshowsPrec = gshowsPrecdefault

instance GShow (f (g p)) => GShow ((f :.: g) p) where
  gshowsPrec = gshowsPrecdefault

instance GShow All where
  gshowsPrec = gshowsPrecdefault

#if MIN_VERSION_base(4,8,0)
instance GShow (f a) => GShow (Alt f a) where
  gshowsPrec = gshowsPrecdefault
#endif

instance GShow Any where
  gshowsPrec = gshowsPrecdefault

#if MIN_VERSION_base(4,9,0)
instance (GShow a, GShow b) => GShow (Arg a b) where
  gshowsPrec = gshowsPrecdefault
#endif

#if !(MIN_VERSION_base(4,9,0))
instance GShow Arity where
  gshowsPrec = gshowsPrecdefault
#endif

instance GShow Associativity where
  gshowsPrec = gshowsPrecdefault

instance GShow Bool where
  gshowsPrec = gshowsPrecdefault

instance GShow BufferMode where
  gshowsPrec = showsPrec

#if defined(HTYPE_CC_T)
instance GShow CCc where
  gshowsPrec = showsPrec
#endif

instance GShow CChar where
  gshowsPrec = showsPrec

instance GShow CClock where
  gshowsPrec = showsPrec

#if defined(HTYPE_DEV_T)
instance GShow CDev where
  gshowsPrec = showsPrec
#endif

instance GShow CDouble where
  gshowsPrec = showsPrec

instance GShow CFloat where
  gshowsPrec = showsPrec

#if defined(HTYPE_GID_T)
instance GShow CGid where
  gshowsPrec = showsPrec
#endif

instance GShow Char where
  gshowsPrec = showsPrec
  gshowList  = showList

#if defined(HTYPE_INO_T)
instance GShow CIno where
  gshowsPrec = showsPrec
#endif

instance GShow CInt where
  gshowsPrec = showsPrec

instance GShow CIntMax where
  gshowsPrec = showsPrec

instance GShow CIntPtr where
  gshowsPrec = showsPrec

instance GShow CLLong where
  gshowsPrec = showsPrec

instance GShow CLong where
  gshowsPrec = showsPrec

#if defined(HTYPE_MODE_T)
instance GShow CMode where
  gshowsPrec = showsPrec
#endif

#if defined(HTYPE_NLINK_T)
instance GShow CNlink where
  gshowsPrec = showsPrec
#endif

#if defined(HTYPE_OFF_T)
instance GShow COff where
  gshowsPrec = showsPrec
#endif

#if MIN_VERSION_base(4,4,0)
instance GShow a => GShow (Complex a) where
  gshowsPrec = gshowsPrecdefault
#endif

instance GShow a => GShow (Const a b) where
  gshowsPrec = gshowsPrecdefault

#if defined(HTYPE_PID_T)
instance GShow CPid where
  gshowsPrec = showsPrec
#endif

instance GShow CPtrdiff where
  gshowsPrec = showsPrec

#if defined(HTYPE_RLIM_T)
instance GShow CRLim where
  gshowsPrec = showsPrec
#endif

instance GShow CSChar where
  gshowsPrec = showsPrec

#if defined(HTYPE_SPEED_T)
instance GShow CSpeed where
  gshowsPrec = showsPrec
#endif

#if MIN_VERSION_base(4,4,0)
instance GShow CSUSeconds where
  gshowsPrec = showsPrec
#endif

instance GShow CShort where
  gshowsPrec = showsPrec

instance GShow CSigAtomic where
  gshowsPrec = showsPrec

instance GShow CSize where
  gshowsPrec = showsPrec

#if defined(HTYPE_SSIZE_T)
instance GShow CSsize where
  gshowsPrec = showsPrec
#endif

#if defined(HTYPE_TCFLAG_T)
instance GShow CTcflag where
  gshowsPrec = showsPrec
#endif

instance GShow CTime where
  gshowsPrec = showsPrec

instance GShow CUChar where
  gshowsPrec = showsPrec

#if defined(HTYPE_UID_T)
instance GShow CUid where
  gshowsPrec = showsPrec
#endif

instance GShow CUInt where
  gshowsPrec = showsPrec

instance GShow CUIntMax where
  gshowsPrec = showsPrec

instance GShow CUIntPtr where
  gshowsPrec = showsPrec

instance GShow CULLong where
  gshowsPrec = showsPrec

instance GShow CULong where
  gshowsPrec = showsPrec

#if MIN_VERSION_base(4,4,0)
instance GShow CUSeconds where
  gshowsPrec = showsPrec
#endif

instance GShow CUShort where
  gshowsPrec = showsPrec

instance GShow CWchar where
  gshowsPrec = showsPrec

instance GShow Double where
  gshowsPrec = showsPrec

instance GShow a => GShow (Down a) where
  gshowsPrec = gshowsPrecdefault

instance GShow a => GShow (Dual a) where
  gshowsPrec = gshowsPrecdefault

instance (GShow a, GShow b) => GShow (Either a b) where
  gshowsPrec = gshowsPrecdefault

instance GShow ExitCode where
  gshowsPrec = gshowsPrecdefault

instance GShow Fd where
  gshowsPrec = showsPrec

instance GShow a => GShow (Monoid.First a) where
  gshowsPrec = gshowsPrecdefault

#if MIN_VERSION_base(4,9,0)
instance GShow a => GShow (Semigroup.First a) where
  gshowsPrec = gshowsPrecdefault
#endif

instance GShow Fixity where
  gshowsPrec = gshowsPrecdefault

instance GShow Float where
  gshowsPrec = showsPrec

instance GShow (ForeignPtr a) where
  gshowsPrec = showsPrec

instance GShow (FunPtr a) where
  gshowsPrec = showsPrec

instance GShow GeneralCategory where
  gshowsPrec = showsPrec

instance GShow Handle where
  gshowsPrec = showsPrec

instance GShow HandlePosn where
  gshowsPrec = showsPrec

#if MIN_VERSION_base(4,8,0)
instance GShow a => GShow (Identity a) where
  gshowsPrec = gshowsPrecdefault
#endif

instance GShow Int where
  gshowsPrec = showsPrec

instance GShow Int8 where
  gshowsPrec = showsPrec

instance GShow Int16 where
  gshowsPrec = showsPrec

instance GShow Int32 where
  gshowsPrec = showsPrec

instance GShow Int64 where
  gshowsPrec = showsPrec

instance GShow Integer where
  gshowsPrec = showsPrec

instance GShow IntPtr where
  gshowsPrec = showsPrec

instance GShow IOError where
  gshowsPrec = showsPrec

instance GShow IOErrorType where
  gshowsPrec = showsPrec

instance GShow IOMode where
  gshowsPrec = showsPrec

instance GShow c => GShow (K1 i c p) where
  gshowsPrec = gshowsPrecdefault

instance GShow a => GShow (Monoid.Last a) where
  gshowsPrec = gshowsPrecdefault

#if MIN_VERSION_base(4,9,0)
instance GShow a => GShow (Semigroup.Last a) where
  gshowsPrec = gshowsPrecdefault
#endif

instance GShow (f p) => GShow (M1 i c f p) where
  gshowsPrec = gshowsPrecdefault

#if MIN_VERSION_base(4,9,0)
instance GShow a => GShow (Max a) where
  gshowsPrec = gshowsPrecdefault
#endif

instance GShow a => GShow (Maybe a) where
  gshowsPrec = gshowsPrecdefault

#if MIN_VERSION_base(4,9,0)
instance GShow a => GShow (Min a) where
  gshowsPrec = gshowsPrecdefault
#endif

#if MIN_VERSION_base(4,8,0)
instance GShow Natural where
  gshowsPrec = showsPrec
#endif

#if MIN_VERSION_base(4,9,0)
instance GShow a => GShow (NonEmpty a) where
  gshowsPrec = gshowsPrecdefault
#endif

instance GShow Ordering where
  gshowsPrec = gshowsPrecdefault

instance GShow p => GShow (Par1 p) where
  gshowsPrec = gshowsPrecdefault

instance GShow a => GShow (Product a) where
  gshowsPrec = gshowsPrecdefault

#if MIN_VERSION_base(4,7,0)
instance GShow (Proxy s) where
  gshowsPrec = gshowsPrecdefault
#endif

instance GShow (Ptr a) where
  gshowsPrec = showsPrec

instance GShow SeekMode where
  gshowsPrec = showsPrec

instance GShow a => GShow (Sum a) where
  gshowsPrec = gshowsPrecdefault

instance GShow (U1 p) where
  gshowsPrec = gshowsPrecdefault

instance GShow (UChar p) where
  gshowsPrec = gshowsPrecdefault

instance GShow (UDouble p) where
  gshowsPrec = gshowsPrecdefault

instance GShow (UFloat p) where
  gshowsPrec = gshowsPrecdefault

instance GShow (UInt p) where
  gshowsPrec = gshowsPrecdefault

instance GShow (UWord p) where
  gshowsPrec = gshowsPrecdefault

instance GShow Version where
  gshowsPrec = gshowsPrecdefault

#if MIN_VERSION_base(4,8,0)
instance GShow Void where
  gshowsPrec = showsPrec
#endif

instance GShow Word where
  gshowsPrec = showsPrec

instance GShow Word8 where
  gshowsPrec = showsPrec

instance GShow Word16 where
  gshowsPrec = showsPrec

instance GShow Word32 where
  gshowsPrec = showsPrec

instance GShow Word64 where
  gshowsPrec = showsPrec

instance GShow WordPtr where
  gshowsPrec = showsPrec

#if MIN_VERSION_base(4,9,0)
instance GShow m => GShow (WrappedMonoid m) where
  gshowsPrec = gshowsPrecdefault
#endif

instance GShow a => GShow (ZipList a) where
  gshowsPrec = gshowsPrecdefault

#if MIN_VERSION_base(4,10,0)
instance GShow CBool where
  gshowsPrec = showsPrec

# if defined(HTYPE_BLKSIZE_T)
instance GShow CBlkSize where
  gshowsPrec = showsPrec
# endif

# if defined(HTYPE_BLKCNT_T)
instance GShow CBlkCnt where
  gshowsPrec = showsPrec
# endif

# if defined(HTYPE_CLOCKID_T)
instance GShow CClockId where
  gshowsPrec = showsPrec
# endif

# if defined(HTYPE_FSBLKCNT_T)
instance GShow CFsBlkCnt where
  gshowsPrec = showsPrec
# endif

# if defined(HTYPE_FSFILCNT_T)
instance GShow CFsFilCnt where
  gshowsPrec = showsPrec
# endif

# if defined(HTYPE_ID_T)
instance GShow CId where
  gshowsPrec = showsPrec
# endif

# if defined(HTYPE_KEY_T)
instance GShow CKey where
  gshowsPrec = showsPrec
# endif

# if defined(HTYPE_TIMER_T)
instance GShow CTimer where
  gshowsPrec = showsPrec
# endif
#endif

max_stars_repo_name: treeowl/linear-generics

max_issues_repo_name: treeowl/linear-generics

max_forks_repo_name: treeowl/linear-generics


text:
module MyPhoto.Actions.Outliers
    ( rmOutliers
    ) where

import           Control.Concurrent ( getNumCapabilities )
import           Control.Concurrent.Async ( mapConcurrently )
import           Control.Concurrent.MSem as MS
import           Control.Monad
import           GHC.Conc ( numCapabilities )
import           System.Console.GetOpt
import           System.IO.Temp
import           System.Directory
import           System.Exit
import           System.FilePath
import           System.Process
import           Graphics.Netpbm
import qualified Data.ByteString as B
import           Data.Maybe (maybe)
import qualified Statistics.Sample as S
import qualified Data.Vector as V (fromList)

import MyPhoto.Model
import MyPhoto.Utils

help :: String
help = let
    header = "Usage: rmoutliers [OPTION...] files..."
  in usageInfo header options

data Options
  = Options
  { optMaxDistance :: Double
  , optSize        :: Int
  , optHelp        :: Bool
  } deriving Show

defaultOptions :: Options
defaultOptions
  = Options
  { optMaxDistance = 200
  , optSize        = 6
  , optHelp        = False
  }

options :: [OptDescr (Options -> Options)]
options =
  [ Option ['d'] ["maxDistance"]
      (OptArg ((\ d opts -> opts { optMaxDistance = d }) . maybe (optMaxDistance defaultOptions) read)
        "DISTANCE")
      "maximum distance to allow"
  , Option ['s'] ["size"]
      (OptArg ((\ s opts -> opts { optSize = s }) . maybe (optSize defaultOptions) read)
        "SIZE")
      "size of reduced image"
  , Option ['h'] ["help"]
      (NoArg (\ opts -> opts { optHelp = True }))
      "print help"
  ]

getMyOpts :: [String] -> IO (Options, [String])
getMyOpts argv = case getOpt Permute options argv of
                   (o,n,[]  ) -> return (foldl (flip id) defaultOptions o, n)
                   (_,_,errs) -> ioError (userError (concat errs ++ help))

computImgVec :: Int -> FilePath -> Img -> IO (Img, [Int])
computImgVec size tmpdir img = do
  let ppmFile = tmpdir </> (takeFileName img ++ ".ppm")

  (_, _, _, pHandle) <- createProcess (proc "convert" [img, "-resize", show size ++ "x" ++ show size ++ "!", ppmFile])
  exitCode <- waitForProcess pHandle
  unless (exitCode == ExitSuccess) $
    fail ("Resize failed with " ++ show exitCode)

  ppmResult <- parsePPM <$> B.readFile ppmFile
  case ppmResult of
    Right ([ppmImg],_) -> let
        imgVec = (pixelDataToIntList . ppmData) ppmImg
      in return (img, imgVec)
    Left err           -> fail ("PPM parsing failed with " ++ err)
    _                  -> fail "PPM parsing failed"

computImgsVecs :: Int -> FilePath -> [Img] -> IO [(Img, [Int])]
computImgsVecs size tmpdir imgs = do
  capabilities <- getNumCapabilities
  sem <- MS.new (capabilities - 2)
  mapConcurrently (MS.with sem . computImgVec size tmpdir) imgs

rmOutliersImpl :: [String] -> [Img] -> PActionBody
rmOutliersImpl _    []            = return (Right [])
rmOutliersImpl args imgs@(img1:_) = let
    dropByDistances :: Double -> [(Img, [Int])] -> [Int] -> IO [(Img, Double)]
    dropByDistances _ []                                 _       = return []
    dropByDistances maxDist ((img, curVec):imgsWithVecs) lastVec = let
        calculateDistance :: [Int] -> [Int] -> Double
        calculateDistance vec1 vec2 =
          sqrt .
          fromIntegral .
          sum .
          map ((\v -> v*v) . uncurry (-)) $ zip vec1 vec2
        dist = calculateDistance lastVec curVec
      in
        if dist < maxDist
        then do
          fmap ((img, dist) :) (dropByDistances maxDist imgsWithVecs curVec)
        else do
          putStrLn ("drop " ++ img ++ " with dist " ++ show dist)
          dropByDistances maxDist imgsWithVecs lastVec
    printStatsOnDists :: [(Img, Double)] -> IO ()
    printStatsOnDists imgsWithoutOutliers = let
        dists = V.fromList (map snd imgsWithoutOutliers)
      in do
        putStrLn $ "Variance = " ++ show (S.variance dists)
        putStrLn $ "Standard deviation = " ++ show (S.stdDev dists)
        putStrLn $ "Mean = " ++ show (S.mean dists)
  in do
    (opts, _) <- getMyOpts args
    print opts
    if optHelp opts
    then return (Left help)
    else withTempDirectory (takeDirectory img1) "_outliers.tmp"
           (\tmpdir -> do
               imgsWithVecs <- computImgsVecs (optSize opts) tmpdir imgs
               imgsWithoutOutliers <- dropByDistances (optMaxDistance opts) imgsWithVecs []
               printStatsOnDists imgsWithoutOutliers
               return (Right (map fst imgsWithoutOutliers))
             )

rmOutliers :: PrePAction
rmOutliers args     = logSeparator "rmoutliers" <> PAction (rmOutliersImpl args)

max_stars_repo_name: maxhbr/myphoto

max_issues_repo_name: maxhbr/myphoto

max_forks_repo_name: maxhbr/myphoto


text:
-- | Helper for writing tests for numeric code
module Test.QuickCheck.Numeric (
    -- * Approximate equality
    eq
  , eqC
    -- * Function monotonicity
  , Monotonicity(..)
  , monotonicFunction
  , monotonicFunctionIEEE
    -- * Inverse function
  , checkInverse
  , checkInverse2
  ) where

import Data.Complex
import qualified Numeric.IEEE    as IEEE

----------------------------------------------------------------
-- Approximate equality
----------------------------------------------------------------

-- | Approximate equality for 'Double'. Doesn't work well for numbers
--   which are almost zero.
eq :: Double                    -- ^ Relative error
   -> Double -> Double -> Bool
eq eps a b
  | a == 0 && b == 0 = True
  | otherwise        = abs (a - b) <= eps * max (abs a) (abs b)

-- | Approximate equality for 'Complex Double'
eqC :: Double                   -- ^ Relative error
    -> Complex Double
    -> Complex Double
    -> Bool
eqC eps a@(ar :+ ai) b@(br :+ bi)
  | a == 0 && b == 0 = True
  | otherwise        = abs (ar - br) <= eps * d
                    && abs (ai - bi) <= eps * d
  where
    d = max (realPart $ abs a) (realPart $ abs b)

----------------------------------------------------------------
-- Function monotonicity
----------------------------------------------------------------

-- | Function monotonicity type.
data Monotonicity
  = StrictInc   -- ^ Strictly increasing function
  | MonotoneInc -- ^ Monotonically increasing function
  | StrictDec   -- ^ Strictly decreasing function
  | MonotoneDec -- ^ Monotonically decreasing function
  deriving (Show,Eq,Ord)

-- | Check that function is nondecreasing. For floating point number
--   it may give spurious failures so 'monotonicFunction'
--   should be used in this case.
monotonicFunction :: (Ord a, Ord b) => Monotonicity -> (a -> b) -> a -> a -> Bool
monotonicFunction cmp f x1 x2
  = f (min x1 x2) `op` f (max x1 x2)
  where
    op = case cmp of
           StrictInc   -> (< )
           MonotoneInc -> (<=)
           StrictDec   -> (> )
           MonotoneDec -> (>=)

-- | Check that function is nondecreasing taking rounding errors into
--   account. This function makes no distinction between strictly
--   increasing function and monotonically increasing function since
--   distinction is pointless for floating point.
--
--   In fact funstion is allowed to decrease less than one ulp in order
--   to guard againist problems with excess precision. On x86 FPU works
--   with 80-bit numbers but doubles are 64-bit so rounding happens
--   whenever values are moved from registers to memory
monotonicFunctionIEEE :: (Ord a, IEEE.IEEE b) => Monotonicity -> (a -> b) -> a -> a -> Bool
monotonicFunctionIEEE cmp f x1 x2
  =  y1 `op` y2
  || abs (y1 - y2) < abs (y2 * IEEE.epsilon)
  where
    y1 = f (min x1 x2)
    y2 = f (max x1 x2)
    op = case cmp of
           StrictInc   -> (<=)
           MonotoneInc -> (<=)
           StrictDec   -> (>=)
           MonotoneDec -> (>=)

----------------------------------------------------------------
-- Function and its inverse
----------------------------------------------------------------

-- | Check that function is inverse. Breaks down near zero.
checkInverse
  :: (Double -> Double) -- ^ Function @f(x)@
  -> (Double -> Double) -- ^ Inverse function @g@, @g(f(x)) = x@
  -> (Double -> Double) -- ^ Derivative of function @f(x)@
  -> Double             -- ^ Relative error for
                        --   @f(x)@. Usually is machine epsilon.
  -> Double             -- ^ Relative error for inverse function
                        --   @g(x)@. Usually is machine epsilon.
  -> Double -> Bool
checkInverse f invF f' eps eps' x
  = x ~= invF y
  where
    (~=) = eq (eps' + abs (y / f' x * eps))
    y    = f x

-- | Check that function is inverse. Breaks down near zero.
checkInverse2
  :: (Double -> Double) -- ^ Function @f(x)@
  -> (Double -> Double) -- ^ Inverse function @g@, @g(f(x)) = x@
  -> (Double -> Double) -- ^ Derivative of function @g(x)@
  -> Double             -- ^ Relative error for
                        --   @f(x)@. Usually is machine epsilon.
  -> Double             -- ^ Relative error for inverse function
                        --   @g(x)@. Usually is machine epsilon.
  -> Double -> Bool
checkInverse2 f invF invF' eps eps' x
  = x ~= invF y
  where
    (~=) = eq (eps' + abs (y * (invF' y * eps)))
    y    = f x

max_stars_repo_name: Shimuuar/quickcheck-numeric

max_issues_repo_name: Shimuuar/quickcheck-numeric

max_forks_repo_name: Shimuuar/quickcheck-numeric


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Graphics.Rendering.Plot.Figure.Plot
-- Copyright   :  (c) A. V. H. McPhail 2010
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- Creation and manipulation of 'Plot's
--
-----------------------------------------------------------------------------

module Graphics.Rendering.Plot.Figure.Plot (
                                            Plot
                                           -- * Plot elements
                                           , Border
                                           , setBorder
                                           , setPlotBackgroundColour
                                           , setPlotPadding
                                           , withHeading
                                           -- * Series data
                                           , D.Abscissa(), D.Ordinate(), D.Dataset()
                                           , SeriesLabel
                                           , D.FormattedSeries()
                                           , D.line, D.point, D.linepoint
                                           , D.impulse, D.step
                                           , D.area
                                           , D.bar
                                           , D.hist
                                           , D.candle, D.whisker
                                           , setDataset
                                           -- * Annotations
                                           , Location, Head, Fill
                                           , AN.arrow
                                           , AN.oval
                                           , AN.rect
                                           , AN.glyph
                                           , AN.text
                                           , AN.cairo
                                           , withAnnotations
                                           -- ** Plot type
                                           , setSeriesType
                                           , setAllSeriesTypes
                                           -- ** Formatting
                                           , D.PlotFormats(..)
                                           , withSeriesFormat
                                           , withAllSeriesFormats
                                           -- * Range
                                           , Scale(..)
                                           , setRange
                                           , setRangeFromData
                                           -- * Axes
                                           , AX.Axis
                                           , AxisType(..),AxisSide(..),AxisPosn(..)
                                           , clearAxes
                                           , clearAxis
                                           , addAxis
                                           , withAxis
                                           -- * BarSetting
                                           , barSetting
                                           -- * SampleData
                                           , sampleData
                                           -- * Legend
                                           , L.Legend
                                           , LegendBorder
                                           , L.LegendLocation(..), L.LegendOrientation(..)
                                           , clearLegend
                                           , setLegend
                                           , withLegendFormat
                                            -- ** Formatting
                                            , Tick(..), TickValues(..), GridLines
                                           , TickFormat(..)
                                           , AX.setTicks
                                           , AX.setGridlines
                                           , AX.setTickLabelFormat
                                           , AX.setTickLabels
                                           , AX.withTickLabelsFormat
                                           , AX.withAxisLabel
                                           , AX.withAxisLine
                                           , AX.withGridLine
                                           ) where

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

--import Data.Eq.Unicode
--import Data.Bool.Unicode
--import Data.Ord.Unicode

import Numeric.LinearAlgebra.Data hiding(Range)

import qualified Data.Array.IArray as A

import Control.Monad.State
import Control.Monad.Reader
--import Control.Monad.Supply

import Prelude hiding(min,max)
import qualified Prelude as Prelude

import Graphics.Rendering.Plot.Types
import Graphics.Rendering.Plot.Defaults
import qualified Graphics.Rendering.Plot.Figure.Text as T
import qualified Graphics.Rendering.Plot.Figure.Plot.Data as D
import qualified Graphics.Rendering.Plot.Figure.Plot.Axis as AX
import qualified Graphics.Rendering.Plot.Figure.Plot.Legend as L
import qualified Graphics.Rendering.Plot.Figure.Plot.Annotation as AN

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

-- | whether to draw a boundary around the plot area
setBorder :: Border -> Plot ()
setBorder b = modify $ \s -> s { _border = b }

-- | set the plot background colour
setPlotBackgroundColour :: Color -> Plot ()
setPlotBackgroundColour c = modify $ \s -> s { _back_colr = c }

-- | set the padding of the subplot
setPlotPadding :: Double -> Double -> Double -> Double -> Plot ()
setPlotPadding l r b t = modify $ \s -> s { _plot_pads = Padding l r b t }

-- | set the heading of the subplot
withHeading :: Text () -> Plot ()
withHeading m = do
                o <- asks _textoptions
                modify $ \s -> s { _heading = execText m o (_heading s) }

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

-- | set the axis range
setRange :: AxisType -> AxisSide -> Scale -> Double -> Double -> Plot ()
setRange XAxis sd sc min max = modify $ \s -> s { _ranges = setXRanges' (_ranges s) }
    where setXRanges' r
              | sc == Log && min <= 0 = error "non-positive logarithmic range"
              | otherwise          = setXRanges sd r
          setXRanges Lower (Ranges (Left _) yr)       = Ranges (Left (Range sc min max)) yr
          setXRanges Lower (Ranges (Right (_,xr)) yr) = Ranges (Right ((Range sc min max,xr))) yr
          setXRanges Upper (Ranges (Left xr) yr)      = Ranges (Right (xr,Range sc min max)) yr
          setXRanges Upper (Ranges (Right (_,xr)) yr) = Ranges (Right (Range sc min max,xr)) yr
setRange YAxis sd sc min max = modify $ \s -> s { _ranges = setYRanges' (_ranges s) }
    where setYRanges' r
              | sc == Log && min <= 0 = error "non-positive logarithmic range"
              | otherwise          = setYRanges sd r
          setYRanges Lower (Ranges xr (Left _))       = Ranges xr (Left (Range sc min max))
          setYRanges Lower (Ranges xr (Right (_,yr))) = Ranges xr (Right ((Range sc min max,yr)))
          setYRanges Upper (Ranges xr (Left yr))      = Ranges xr (Right (yr,Range sc min max))
          setYRanges Upper (Ranges xr (Right (_,yr))) = Ranges xr (Right ((Range sc min max,yr)))

-- | set the axis ranges to values based on dataset
setRangeFromData :: AxisType -> AxisSide -> Scale -> Plot ()
setRangeFromData ax sd sc = do
  ds <- gets _data
  let ((xmin,xmax),(ymin,ymax)) = calculateRanges ds
  case ax of
    XAxis -> setRange ax sd sc (if sc == Log then if xmin == 0 then 1 else xmin else xmin) xmax
    YAxis -> setRange ax sd sc (if sc == Log then if ymin == 0 then 1 else ymin else ymin) ymax

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

withAnnotations :: Annote () -> Plot ()
withAnnotations = annoteInPlot

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

-- | clear the axes of a subplot
clearAxes :: Plot ()
clearAxes = modify $ \s -> s { _axes = [] }

-- | clear an axis of a subplot
clearAxis :: AxisType -> AxisPosn -> Plot ()
clearAxis at axp = do
  ax <- gets _axes
  modify $ \s -> s { _axes =
    filter (\(Axis at' axp' _ _ _ _ _ _) -> not (at == at' && axp == axp')) ax }

-- | add an axis to the subplot
addAxis :: AxisType -> AxisPosn -> AX.Axis () -> Plot ()
addAxis at axp m  = do
  ax' <- gets _axes
  o <- ask
  let ax = execAxis m o (defaultAxis at axp)
  modify $ \s -> s { _axes = ax : ax' }

-- | operate on the given axis
withAxis :: AxisType -> AxisPosn -> AX.Axis () -> Plot ()
withAxis at axp m = do
  axes' <- gets _axes
  o <- ask
  modify $ \s -> s { _axes =
    map (\a@(Axis at' ap' _ _ _ _ _ _) -> if at == at' && axp == ap'
                                         then execAxis m o a
                                         else a) axes' }

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

barSetting :: BarSetting -> Plot ()
barSetting bc = modify $ \s -> s { _barconfig = bc }

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

sampleData :: SampleData -> Plot ()
sampleData sd = modify $ \s -> s { _sampledata = sd }

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

-- | clear the legend
clearLegend :: Plot ()
clearLegend = withLegend $ L.clearLegend

-- | set the legend location and orientation
setLegend :: L.LegendBorder -> L.LegendLocation -> L.LegendOrientation -> Plot ()
setLegend b l o = withLegend $ L.setLegend b l o

-- | format the legend text
withLegendFormat :: T.Text () -> Plot ()
withLegendFormat f = withLegend $ L.withLegendFormat f

-- | operate on the legend
withLegend :: L.Legend () -> Plot ()
withLegend = legendInPlot

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

-- | operate on the data
withData :: D.Data () -> Plot ()
withData = dataInPlot

{- | set the data series of the subplot

   The data series are either 'FormattedSeries' or plain data series.
   A plain data series must carry a 'SeriesType'.

  A dataset may or may not have an abscissa series, and if so, it is paired
  with either a list of ordinate series or a single ordinate series.

  The abscissa series (if present) is of type 'Vector Double'.

  An ordinate series be a function (@Double -> Double@) or a series of points,
  a 'Vector Double' with optional error series, y axis preference, and labels.

  To specify decoration options for an ordinate series, use the appropriate function, such
  as 'linespoints', with the ordinate series and decoration formatting ('LineFormat',
  'PointFormat', and 'BarFormat') as arguments.

> setDataset (ts,[linespoints (xs,(le,ue),Upper,"data") (([Dash,Dash],3,blue),(Diamond,green))])

  has abscissa @ts@ paired with a list of ordinate series, the single element of which is a
  'FormattedSeries', @linespoints@ where the ordinate is @xs@ with error series @le@ and @ue@,
  to be graphed against the upper y-range with label \"data\".  The line element is formatted
  to be dashed, of width 3, and blue and the point element is to be a green diamond.
-}
setDataset :: D.Dataset a => a -> Plot ()
setDataset d = withData $ D.setDataSeries d

-- | set the plot type of a given data series
setSeriesType :: Int -> SeriesType -> Plot ()
setSeriesType i t = withData $ D.setSeriesType t i

-- | change the plot type of all data series
setAllSeriesTypes :: SeriesType -> Plot ()
setAllSeriesTypes t = withData $ D.setAllSeriesTypes t

-- | format the plot elements of a given series
withSeriesFormat :: D.PlotFormats m => Int -> m () -> Plot ()
withSeriesFormat i f = withData $ D.withSeriesFormat i f

{- |
  format the plot elements of all series

      the operation to modify the formats is passed the series index.
      This allows, for example, colours to be selected from a list
      that gets indexed by the argument

> setColour = withAllSeriesFormats (\i -> do
>                                         setLineColour $ [black,blue,red,green,yellow] !! i
>                                         setLineWidth 1.0)
-}
withAllSeriesFormats :: D.PlotFormats m => (Int -> m ()) -> Plot ()
withAllSeriesFormats f = withData $ D.withAllSeriesFormats f

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

findMinMax :: Abscissae -> Ordinates -> (Double,Double)
findMinMax (AbsFunction _) (OrdFunction _ f _) =
    let v = f (linspace 100 (-1,1))
    in (minElement v,maxElement v)
findMinMax (AbsPoints _ x) (OrdFunction _ f _) =
    let v = f x
    in (minElement v,maxElement v)
-- what if errors go beyond plot?
findMinMax _ (OrdPoints _ (Plain o) _)    = (minElement o,maxElement o)
findMinMax _ (OrdPoints _ (Error o _) _)  = (minElement o,maxElement o)
findMinMax _ (OrdPoints _ (MinMax (o,p) _) _) =
  (Prelude.min (minElement o) (minElement p)
  ,Prelude.max (maxElement o) (maxElement p))

abscMinMax :: Abscissae -> (Double,Double)
abscMinMax (AbsFunction _) = defaultXAxisSideLowerRange
abscMinMax (AbsPoints _ x) = (minElement x,maxElement x)

ordDim :: Ordinates -> Int
ordDim (OrdFunction _ _ _)  = 1
ordDim (OrdPoints _ o _)    = size $ getOrdData o

calculateRanges :: DataSeries -> ((Double,Double),(Double,Double))
calculateRanges (DS_Y ys)      =
  let xmax = maximum $ map (\(DecSeries o _) ->
        fromIntegral $ ordDim o) $ A.elems ys
      ym = unzip $ map (\(DecSeries o _) ->
        findMinMax (AbsFunction id) o) $ A.elems ys
      ymm = (minimum $ fst ym,maximum $ snd ym)
  in ((0,xmax),ymm)
calculateRanges (DS_1toN x ys) =
  let ym = unzip $ map (\(DecSeries o _) -> findMinMax x o) $ A.elems ys
      ymm = (minimum $ fst ym,maximum $ snd ym)
      xmm = abscMinMax x
  in (xmm,ymm)
calculateRanges (DS_1to1 ys)   =
  let (xm',ym') = unzip $ A.elems ys
      ym = unzip $ map (\(x,(DecSeries o _)) -> findMinMax x o) (zip xm' ym')
      ymm = (minimum $ fst ym,maximum $ snd ym)
      xm = unzip $ map abscMinMax xm'
      xmm = (minimum $ fst xm,maximum $ snd xm)
  in (xmm,ymm)
calculateRanges (DS_Surf m)     =
    ((0,fromIntegral $ cols m),(fromIntegral $ rows m,0))

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

max_stars_repo_name: JackTheEngineer/plot


text:
module Space where

import Data.Set
import qualified Data.Set as S

import Data.Vector
import qualified Data.Vector as V

import Numeric.LinearAlgebra
import qualified Numeric.LinearAlgebra as M

-- INITIALISATION ---------------------------------------------------------------------------------

-- The data Space (short for simplicial complex) consist of a set of simplices with
-- the following property: all the faces of all the simplices in the set must also
-- belong to the set. A simplex is implemented as a Set for efficiency.
-- For convenience, the set of simplices is split into a list according to the
-- dimension. The defining condition is unchecked as, in practice, the helper
-- function allSimplices (below) is always used to iniatialise the data. The output
-- of the helper function is always valid.

data Space t = Space { subsimplices :: [Set (Set t)] } deriving (Show)

-- This helper function takes a set of simplices as input and returns a valid
-- simplicial complex. The input consists of all the highest dimensional simplices
-- which make up the space. The function takes all their boundaries and sorts
-- them by dimension. The lowest dimension is -1 and it is the dimension of the
-- empty set. It is considered a simplex.
allSimplices :: (Ord t) => [[t]] -> [Set (Set t)]
allSimplices [] = error "This space has no faces. For empty space input '[[]]'."
allSimplices faces =

    let associatedSet = S.fromList (Prelude.map (S.fromList) faces)
        powerFaces = S.map (S.powerSet) associatedSet
        allMixedFaces = S.foldl (S.union) S.empty powerFaces
        lengthFunctions = [ \xs -> (S.size xs == n) | n <- [0,1..] ]

    in  sieveAndForget lengthFunctions allMixedFaces

sieveAndForget :: [t -> Bool] -> Set t -> [Set t]
sieveAndForget (f:fs) xs
    | (S.size xs == 0) = []
    | otherwise        = p1 : (sieveAndForget fs p2)
        where (p1,p2) = S.partition (f) xs

-- HOMOLOGY ---------------------------------------------------------------------------------------

-- This function only computes the free rank of the homology groups. It does not give
-- any information about the torsion subgroups. These require knowledge of the
-- Smith decomposition of a matrix, an algorithm which has never been developed in
-- Haskell (over a PID). To compute the free rank, it is enough to know the rank
-- of the boundary operators (through the package Numeric.LinearALgebra)

-- This function returns the homology numbers of a space from 0 up to its dimension.
-- All the other homology numbers must be zero.
homology :: (Ord t) => Space t -> [Int]
homology (Space { subsimplices = xs }) =
    let operators = chainComplex (Space { subsimplices = xs })
        dimensionDomains = Prelude.map (M.cols) operators
        dimensionImages = Prelude.map (M.rank) (Prelude.tail operators)
        dimensionAugmentedImages = 0:dimensionImages
        dimensionShiftedImages = dimensionImages Prelude.++ [0]
        dimensionKernels = Prelude.zipWith (-) dimensionDomains dimensionAugmentedImages

    in  Prelude.zipWith (-) dimensionKernels dimensionShiftedImages

-- Very small variation on homology. It is used to understand if a space is connected.
reducedHomology :: (Ord t) => Space t -> [Int]
reducedHomology space = reduced
    where reduced = (x-1):xs
          (x:xs) = homology space

-- Computes all boundary operators in a chain complex
chainComplex :: (Ord t) => Space t -> [M.Matrix R]
chainComplex (Space { subsimplices = xs }) =

    let listOfSimplices = Prelude.map (S.toList) xs
        domainListOfSimplices = Prelude.tail listOfSimplices
        imageListOfSimplices = Prelude.init listOfSimplices
        domainVectorOfSimplices = Prelude.map (V.fromList) domainListOfSimplices
        imageVectorOfSimplices = Prelude.map (V.fromList) imageListOfSimplices
        domainImage = Prelude.zip domainVectorOfSimplices imageVectorOfSimplices

    in  Prelude.map (boundaryOperator) domainImage

-- Computes boundary operator for all simplices of a fixed dimension.
boundaryOperator :: (Ord t) => (V.Vector (Set t) , V.Vector (Set t))-> M.Matrix R
boundaryOperator (domain , image) = convertToMatrix vectorMatrix
    where vectorMatrix = V.map (f) domain
          f = \x -> (piecewiseBoundaryOperator x image)

-- Computes boundary operator of a single simplex.
piecewiseBoundaryOperator :: (Ord t) => Set t -> V.Vector (Set t) -> V.Vector Int
piecewiseBoundaryOperator element image =
    let f imageFace
              | (what == False) = 0
              | odd indicator   = 1
              | otherwise       = -1

              where (s1 , what , s2) = S.splitMember imageFace facesOfElement
                    -- Shortcut to compute the sign of a face
                    indicator = S.size element + S.size s1
                    facesOfElement = S.map (g) element
                    g = \y -> (S.delete y element)

    in  V.map (f) image

-- MATRIX CONVERSION ------------------------------------------------------------------------------

-- These algorithms are necessary to convert from V.Vector to M.Matrix.
-- Unfortunately these two data lack compatibility. Furthermore
-- M.Vector cannot store Sets, hence they cannot be used from the beginning.
-- Ideally, the Smith Normal Form should be called directly on the V.Vector
-- so to avoid using M.Matrix all together.

convertToMatrix :: V.Vector (V.Vector Int) -> M.Matrix R
convertToMatrix vectorMatrix = M.matrix s doubleFlattenedListMatrix
    where s = V.length vectorMatrix
          doubleFlattenedListMatrix = Prelude.map (fromIntegral) flattenedListMatrix
          flattenedListMatrix = Prelude.foldl (Prelude.++) [] listMatrix
          listMatrix = Prelude.map (V.toList) listVectorMatrix
          listVectorMatrix = V.toList transposedVectorMatrix
          transposedVectorMatrix = transpose vectorMatrix

transpose :: V.Vector (V.Vector Int) -> V.Vector (V.Vector Int)
transpose matrix = V.foldr (g) emptyVector matrix
    where emptyVector = V.replicate s V.empty
          s = V.length ((V.!) matrix 0)
          g = \x y -> V.zipWith (V.cons) x y

-- BASIC FUNCTIONS --------------------------------------------------------------------------------

-- These functions are basic as they do not rely on computing homology.

-- The function returns the dimension of the highest dimensional simplices.
-- 2 is subtracted to account for the simplices of dimension 0 and -1.
dimension :: (Ord t) => Space t -> Int
dimension (Space { subsimplices = xs }) = Prelude.length xs - 2

-- This function returns the number of simplices of a given dimension.
numberOfSimplices :: (Ord t) => Space t -> Int -> Int
numberOfSimplices (Space { subsimplices = xs }) dim
    | dim < -1 || dim > maxDim = 0
    | otherwise                   = S.size simplicesDim

    where maxDim = dimension (Space { subsimplices = xs })
          simplicesDim = xs !! index
          index = dim + 1

eulerCharacteristic :: (Ord t) => Space t -> Int
eulerCharacteristic (Space { subsimplices = xs }) =
    alternatingSum (Prelude.tail lengths)

        where alternatingSum :: (Num t) => [t] -> t
              alternatingSum [] = error "Not defined."
              alternatingSum [x] = x
              alternatingSum (x:xs) = x - alternatingSum xs

              lengths = (Prelude.map (S.size) xs)

-- ADVANCED FUNCTIONS -----------------------------------------------------------------------------

-- These functions all rely on homology.

isConnected :: (Ord t) => Space t -> Bool
isConnected space = (zeroHomology == 1)
    where zeroHomology = Prelude.head (homology space)

-- This function can tell if a space is not homotopic to a topological manifold.
-- Note that, for the reasons above, it only checks Poincaré duality up to torsion.
rationalPoincareDuality :: (Ord t) => Space t -> Bool
rationalPoincareDuality space = (homNumbers == revHomNumbers)
    where revHomNumbers = Prelude.reverse homNumbers
          homNumbers = homology space

-- If this function returns False then the answer is definitely No. If it returns
-- True, then the algorthm is inconcludent. The function first checks the Euler
-- characteristic as this is a less expensice invariant to compute, albeit coarser.
-- If it is the same, homology ranks are computed.
couldItBeContractible :: (Ord t) => Space t -> Bool
couldItBeContractible space = (euler == 1) && (isZero reduced)
    where euler = eulerCharacteristic space
          reduced = reducedHomology space

-- Simple helper function to check if all the entries in a list are zero
isZero :: [Int] -> Bool
isZero [] = error "Empty list."
isZero [0] = True
isZero (x:xs) = (x == 0) && (isZero xs)

-- Again, if the function returns True, no new information is gained. The
-- Euler characteristic is computed first. This algorithm only checks
-- if the homology numbers are the same, however, it is necessary to add
-- enough zeros to the shorter list to equate the lengths.
couldTheyBeHomotopic :: (Ord t) => Space t -> Space t -> Bool
couldTheyBeHomotopic space1 space2 = (euler1 == euler2) && (hom1 == hom2)
    where euler1 = eulerCharacteristic space1
          euler2 = eulerCharacteristic space2
          hom1 = tempHom1 Prelude.++ zeroList1
          hom2 = tempHom2 Prelude.++ zeroList2
          tempHom1 = homology space1
          tempHom2 = homology space2
          zeroList1 = [ 0 | x <- [1..s1]]
          zeroList2 = [ 0 | x <- [1..s2]]
          s1 = maxLength - length1
          s2 = maxLength - length2
          length1 = Prelude.length tempHom1
          length2 = Prelude.length tempHom2
          maxLength = Prelude.maximum [length1, length2]


text:
{-# LANGUAGE LambdaCase #-}
module Lib
( main
) where

import Prelude hiding (lex)
import Control.Applicative (Alternative((<|>)))
import Numeric.LinearAlgebra
    (realPart, Complex(..), eigenvalues, (><), toList, Matrix)
import Data.Attoparsec.Text
    (many', endOfInput, double, skipSpace, parseOnly, Parser)
import Data.Text (pack)

main :: IO ()
main = do
  input <- fmap (\case {';' -> ' '; ch -> ch}) <$> getContents
  let result = parseOnly (matrixParser <* endOfInput) $ pack input
  putStrLn $ case result of
    Left  err -> err
    Right m   -> definiteness m

matrixParser :: Parser (Matrix (Complex Double))
matrixParser = loadMatrix =<< (skipSpace *> many' (double <* skipSpace))
  where
    loadMatrix xs = do
      let n = floor . sqrt . fromIntegral . length $ xs
      if n * n /= length xs
      then fail "The provided matrix has to be square"
      else pure $ (n><n) $ (:+ 0) <$> xs

definiteness :: Matrix (Complex Double) -> String
definiteness m = flip orElse "indefinite" $
  foldl ((. uncurry classify) <$> (<|>)) Nothing classes

  where
    eValues = realPart <$> toList (eigenvalues m)
    classify f c = if all f eValues then Just c else Nothing

    orElse Nothing  x = x
    orElse (Just x) _ = x

    classes = [ ((>  0), "positive definite")
              , ((<  0), "negative definite")
              , ((>= 0), "positive semidefinite")
              , ((<= 0), "negative semidefinite")
              ]

max_stars_repo_name: viluon/mpi-definite-matrices

max_issues_repo_name: viluon/mpi-definite-matrices

max_forks_repo_name: viluon/mpi-definite-matrices


text:
module Statistics.GModeling.Models.HMM
  (
  ) where

import Statistics.GModeling.DSL
import Statistics.GModeling.Gibbs
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed.Mutable as MU

data HMM = Alpha | Beta | Transition
         | Topic | Symbols | Symbol

hmm :: Network HMM
hmm =
  [
    Only Alpha            :-> Transition
  , Only Beta             :-> Symbols
  , (Transition :@ Topic) :-> Topic
  , (Symbols :@ Topic)    :-> Symbol
  ]

type D = V.Vector (U.Vector (Int,Int))

reader :: D -> Reader HMM Int
reader v = Reader
  {
    size = G.sum (G.map G.length v)
  , readn = \idx k ->
      let (i,j) = indices G.! idx
          (topic,symbol) = v G.! i G.! j
          (prev_topic,_) = v G.! i G.! (j-1)
      in case k of
           Topic -> topic
           Symbol -> symbol
           Symbols -> topic
           Transition -> if j==0
                         -- reserve Topic 0 for initial distribution
                         then 0
                         else prev_topic
  , copy = do
      m <- G.mapM U.thaw v
      return Writer
        {
          writen = \idx k v -> do
            let (i,j) = indices G.! idx
            case k of
              Topic -> do
                (_,s) <- MU.unsafeRead (m G.! i) j
                MU.unsafeWrite (m G.! i) j (v,s)
        , readOnly = G.mapM U.freeze m >>= return . reader
        }
  }
  where indices = U.fromList $ do
          i <- [0..G.length v-1]
          j <- [0..G.length (v G.! i)-1]
          return (i,j)

max_stars_repo_name: nanonaren/gmodeling

max_issues_repo_name: nanonaren/gmodeling

max_forks_repo_name: nanonaren/gmodeling


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

module Utils where
import Control.Monad.Fix
import Control.Monad
import System.IO
import Signal
import Circuit
import AmpOp
    ( SignalState(SignalTwoNodes),
      AmpOp(openLoopGain),
      Impedance,
      Output,
      Input,
      dc100,
      sen,
      ground,
      lm741,
      r1,
      r2,
      eps )
import Data.List
import Data.Complex
import Simulation

type Iteration = Int

makePlotIterations = do
    outh <- openFile "outputIterations.txt" WriteMode
    mapM_ (\(i, o) -> hPutStrLn outh (show i ++ " " ++ show o)) (ampOpInvertingList' lm741 12 r1 r2)
    hClose outh

ampOpInvertingList' :: AmpOp -> Input -> Impedance -> Impedance -> [(Iteration, Output)]
ampOpInvertingList' model input r1 r2 = fix calculate (100, 0, [], 1)
    where calculate f (vXOld, vOutOld, l, it) = if abs(realPart vOutOld - realPart vOut) <= 0.000001 && abs(realPart vXOld - realPart vX) <= 0.000001 then l ++ [(it, vOut)] else f (vOut, vX, l ++ [(it, vOut)], it + 1)
            where vX = (- vOutOld) / openLoopGain model `at` 0
                  vOut = (vXOld * (r1 + r2) `at` 0 - input * r2 `at` 0) / r1 `at` 0

checkDataIterations :: IO ()
checkDataIterations = print (ampOpInvertingList' lm741 12 r1 r2)

ampOpInvertingList :: AmpOp -> Impedance -> Impedance -> SignalState -> Circuit Input [(Iteration, SignalState)]
ampOpInvertingList model r1 r2 (SignalTwoNodes (vXInit, vOutInit)) =
    ($ (vXInit, vOutInit, [], 1))<$> mfix
    (\f -> Circuit $ \vIn -> Signal $ \time (vXOld, vOutOld, l, it) ->
        let vX = (- vOutOld) / openLoopGain model
            vOut = ((vXOld * (r1 + r2)) - (vIn * r2)) / r1
        in if abs((realPart <$> vOutOld) - (realPart <$> vOut)) <= eps && abs((realPart <$> vXOld) - (realPart <$> vX)) <= eps
            then l ++ [(it, SignalTwoNodes (vX, vOut))]
            else f (vX, vOut, l ++ [(it, SignalTwoNodes (vX, vOut))], it + 1))

makeSamplesInput :: Signal Input -> [Metric] -> [(Metric, Input)]
makeSamplesInput input samples = do time <- samples
                                    let result = input `at` time
                                    return (time, result)

makeSamplesOutput :: SignalState -> (SignalState -> Circuit Input [(Iteration, SignalState)]) -> Signal Input -> [Metric] -> [[(Metric, Iteration, Output)]]
makeSamplesOutput _ _ _ [] = []
makeSamplesOutput initialState circuit input (t:ts) = map (\(i, SignalTwoNodes (_, o)) -> (t, i, o `at` t)) results : rest
    where results = circuit initialState `simulate` input `at` t
          rest = makeSamplesOutput (snd $ last results) circuit input ts

ampOpInvertingList2 :: AmpOp -> Impedance -> Impedance -> Circuit Input [(Iteration, Output)]
ampOpInvertingList2 model r1 r2 =
    ($ (dc100, ground, [], 1))<$> mfix
    (\f -> Circuit $ \vIn -> Signal $ \time (vXOld, vOutOld, l, it) ->
        let vX = (- vOutOld) / openLoopGain model
            vOut = ((vXOld * (r1 + r2)) - (vIn * r2)) / r1
        in if abs((realPart <$> vOutOld) - (realPart <$> vOut)) <= eps && abs((realPart <$> vXOld) - (realPart <$> vX)) <= eps
            then l ++ [(it, vOut `at` time)]
            else f (vX, vOut, l ++ [(it, vOut `at` time)], it + 1))

makeSamplesOutput2 :: Signal Input -> [Metric] -> [[(Metric, Iteration, Output)]]
makeSamplesOutput2 input samples = do time <- samples
                                      let results = (ampOpInvertingList2 lm741 r1 r2 `simulate` input) `at` time
                                      mapM (\(i, o) -> [(time, i, o)]) results

testCircuit = ampOpInvertingList lm741 r1 r2

makePlotMetricOutput = do
    outh <- openFile "outputMetric.txt" WriteMode
    mapM_ (\(t, i, o) -> hPutStrLn outh (show t ++ " " ++ show i ++ " " ++ show o)) (concat $ makeSamplesOutput iSignalState testCircuit sen (map (:+ 0) [0.1,0.2..10.0]))
    hClose outh

makePlotMetricInput = do
    outh <- openFile "inputMetric.txt" WriteMode
    mapM_ (\(t, o) -> hPutStrLn outh (show t ++ " " ++ show o)) (makeSamplesInput sen (map (:+ 0) [0.1,0.2..10.0]))
    hClose outh

max_stars_repo_name: FP-Modeling/fixingAnalog

max_issues_repo_name: FP-Modeling/fixingAnalog

max_forks_repo_name: FP-Modeling/fixingAnalog


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

module Main where

import Data.Array
import Data.Complex
import qualified Data.IntMap.Strict as IM
import Data.List (intercalate)
import Data.Map (empty, fromList, union)
import Data.Maybe (fromJust)
import Data.STRef.Strict
import qualified Data.Set as Set
import HashedExpression.Derivative.Partial
import Graphics.EasyPlot
import HashedExpression
import HashedExpression.Derivative
import HashedExpression.Interp
import HashedExpression.Operation
import qualified HashedExpression.Operation
import HashedExpression.Prettify
import Data.String.Interpolate
import Prelude hiding ((^))

main :: IO ()
main = do
  let x = variable "x"
  let y = variable "y"
  let f = sum [x^2 , 2*x , 1] * (x^2 + 2 * y + 1)
  print $ allEntries f

max_stars_repo_name: Turboscient/HashedExpression

max_issues_repo_name: Turboscient/HashedExpression

max_forks_repo_name: Turboscient/HashedExpression


text:
{-# LANGUAGE BangPatterns,RecordWildCards,FlexibleContexts,TypeFamilies #-}
{-# OPTIONS_GHC -Wno-deferred-type-errors #-}

module Neuro where

{-@ LIQUID "--no-termination" @-}
{-@ LIQUID "--reflection" @-}
{-@ LIQUID "--no-total" @-}

import Codec.Compression.Zlib     (compress, decompress)

import Data.Binary                (Binary(..), encode, decode)
import Data.List                  (foldl')
import Foreign.Storable           (Storable)

import qualified Data.ByteString.Lazy  as BS
import qualified Data.Vector           as DV
      (Vector,
      toList,
      fromList,
      snoc,empty,
      foldl',
      scanl,scanr,
      zipWith,
      tail,init,last,zip,map, head)
import Numeric.LinearAlgebra
    ( (#>),
      reshape,
      vector,
      cmap,
      (<.>),
      (><),
      outer,
      vjoin,
      Element,
      Matrix,
      Container,
      Linear(scale),
      Numeric,
      RealOf,
      Transposable(tr),
      Normed(norm_2),
      Vector)
import System.Random.MWC as MWC
    ( Variate, asGenST, uniformVector, withSystemRandom )

newtype Network a = Network
                 { matrices   :: DV.Vector (Matrix a)
                 } deriving Show

instance (Element a, Binary a) => Binary (Network a) where
  put (Network ms) = put . DV.toList $ ms
  get = (Network . DV.fromList) `fmap` get

type ActivationFunction a = a -> a

type ActivationFunctionDerivative a = a -> a

trainNTimes :: (Floating (Vector a), Floating a, Numeric a, Num (Vector a), Container Vector a) => Int -> a -> ActivationFunction a -> ActivationFunctionDerivative a -> Network a -> Samples a -> Network a
trainNTimes n = trainUntil (\k _ _ -> k > n)

trainUntil :: (Floating (Vector a), Floating a, Numeric a, Num (Vector a), Container Vector a) => (Int -> Network a -> Samples a -> Bool) -> a -> ActivationFunction a -> ActivationFunctionDerivative a -> Network a -> Samples a -> Network a
trainUntil pr learningRate act act' net samples = go net 0
  where go n !k | pr k n samples = n
                | otherwise      = case backpropOnce learningRate act act' n samples of
                                    n' -> go n' (k+1)

createNetwork :: (Variate a, Storable a) => Int -> [Int] -> Int -> IO (Network a)
createNetwork nInputs hiddens nOutputs =
  fmap Network $ withSystemRandom . asGenST $ \gen -> go gen dimensions DV.empty
  where
        go _ [] !ms         = return ms
        go gen ((!n,!m):ds) ms = do
          !mat <- randomMat n m gen
          go gen ds (ms `DV.snoc` mat)
        randomMat n m g = reshape m `fmap` uniformVector g (n*m)
        dimensions      = zip (hiddens ++ [nOutputs]) $
                              (nInputs+1 : hiddens)

fromWeightMatrices :: Storable a => DV.Vector (Matrix a) -> Network a
fromWeightMatrices ws = Network ws

output :: (Floating (Vector a), Numeric a, Storable a, Num (Vector a)) => Network a -> ActivationFunction a -> Vector a -> Vector a
output (Network{..}) act input = DV.foldl' f (vjoin [input, 1]) matrices
  where f !inp m = cmap act $ m #> inp

outputs :: (Floating (Vector a), Numeric a, Storable a, Num (Vector a)) => Network a -> ActivationFunction a -> Vector a -> DV.Vector (Vector a)
outputs (Network{..}) act input = DV.scanl f (vjoin [input, 1]) matrices
  where f !inp m = cmap act $ m #> inp

deltas :: (Floating (Vector a), Floating a, Numeric a, Container Vector a, Num (Vector a)) => Network a -> ActivationFunctionDerivative a -> DV.Vector (Vector a) -> Vector a -> DV.Vector (Matrix a)
deltas (Network{..}) act' os expected = DV.zipWith outer (DV.tail ds) (DV.init os)
  where !dl = (DV.last os - expected) * (deriv $ DV.last os)
        !ds = DV.scanr f dl (DV.zip os matrices)
        f (!o, m) !del = deriv o * (tr m #> del)
        deriv = cmap act'

updateNetwork :: (Floating (Vector a), Floating a, Numeric a, Storable a, Num (Vector a), Container Vector a) => a -> ActivationFunction a -> ActivationFunctionDerivative a -> Network a -> Sample a -> Network a
updateNetwork alpha act act' n@(Network{..}) (input, expectedOutput) = Network $ DV.zipWith (+) matrices corr
    where !xs = outputs n act input
          !ds = deltas n act' xs expectedOutput
          !corr = DV.map (scale (-alpha)) ds

type Sample a = (Vector a, Vector a)

type Samples a = [Sample a]

backpropOnce :: (Floating (Vector a), Floating a, Numeric a, Num (Vector a), Container Vector a) => a -> ActivationFunction a -> ActivationFunctionDerivative a -> Network a -> Samples a -> Network a
backpropOnce rate act act' n samples = foldl' (updateNetwork rate act act') n samples

quadError :: (Floating (Vector a), Floating a, Fractional (RealOf a), Normed (Vector a), Numeric a) => ActivationFunction a -> Network a -> Samples a -> RealOf a
quadError act net samples = realToFrac $ foldl' (\err (inp, out) -> err + (norm_2 $ output net act inp - out)) 0 samples

tanh' :: Floating a => a -> a
tanh' x = let s = tanh x
           in 1 - s**2

loadNetwork :: (Storable a, Element a, Binary a) => FilePath -> IO (Network a)
loadNetwork fp = return . decode . decompress =<< BS.readFile fp

saveNetwork :: (Storable a, Element a, Binary a) => FilePath -> Network a -> IO ()
saveNetwork fp net = BS.writeFile fp . compress $ encode net


text:
module STCPinwheelCorner where

import           Control.Monad                  as M
import           Control.Monad.Parallel         as MP
import qualified Data.Array.Accelerate          as A
import           Data.Array.Accelerate.LLVM.PTX
import           Data.Array.IArray              as IA
import           Data.Array.Repa                as R
import           Data.Binary
import           Data.Complex
import           Data.List                      as L
import           Data.Vector.Storable           as VS
import           Data.Vector.Unboxed            as VU
import           Filter.Utils
import           FokkerPlanck
import           Foreign.CUDA.Driver            as CUDA
import           FourierMethod.BlockMatrixAcc
import           FourierMethod.FourierSeries2D
import           FourierPinwheel
import           Image.IO
import           Pinwheel.FourierSeries2D
import           STC                            hiding (convolve)
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.Array
import           Utils.Time

main = do
  args@(deviceIDsStr:numPointsStr:deltaStr:thresholdStr:numPointsReconStr:deltaReconStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:tauStr:numR2FreqStr:periodR2Str:phiFreqsStr:rhoFreqsStr:thetaFreqsStr:scaleFreqsStr:initDistStr:initScaleStr:histFilePath:stdR2Str:stdStr:numBatchR2Str:numBatchR2FreqsStr:numBatchOriStr:batchSizeStr:sStr:radiusStr:deltaTStr:weightStr:posStr:numThreadStr:_) <-
    getArgs
  let deviceIDs = read deviceIDsStr :: [Int]
      numPoints = read numPointsStr :: Int
      delta = read deltaStr :: Double
      threshold = read thresholdStr :: Double
      numPointsRecon = read numPointsReconStr :: Int
      deltaRecon = read deltaReconStr :: Double
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      tau = read tauStr :: Double
      numR2Freq = read numR2FreqStr :: Int
      periodR2 = read periodR2Str :: Double
      phiFreq = read phiFreqsStr :: Int
      phiFreqs = L.map fromIntegral [-phiFreq .. phiFreq]
      rhoFreq = read rhoFreqsStr :: Int
      rhoFreqs = L.map fromIntegral [-rhoFreq .. rhoFreq]
      thetaFreq = read thetaFreqsStr :: Int
      thetaFreqs = L.map fromIntegral [-thetaFreq .. thetaFreq]
      scaleFreq = read scaleFreqsStr :: Int
      scaleFreqs = L.map fromIntegral [-scaleFreq .. scaleFreq]
      initScale = read initScaleStr :: Double
      initDist = read initDistStr :: [(Double, Double, Double, Double)]
      initPoints = L.map (\(x, y, t, s) -> Point x y t s) initDist
      initSource = [L.head initPoints]
      initSink = [L.last initPoints]
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCPinwheelCorner"
      stdR2 = read stdR2Str :: Double
      std = read stdStr :: Double
      numBatchR2 = read numBatchR2Str :: Int
      numBatchR2Freqs = read numBatchR2FreqsStr :: Int
      numBatchOri = read numBatchOriStr :: Int
      batchSize = read batchSizeStr :: Int
      s = read sStr :: Double
      radius = read radiusStr :: Double
      deltaT = read deltaTStr :: Double
      weight = read weightStr :: Double
      pos = read posStr :: (Int, Int)
      periodEnv = periodR2 ^ 2 / 4
  -- removePathForcibly folderPath
  createDirectoryIfMissing True folderPath
  -- let deltaTheta = 2 * pi / fromIntegral numOrientation
  --     cornerDist =
  --       cornerDistribution
  --         delta
  --         initScale
  --         thetaSigma
  --         tau
  --         threshold
  --         [fromIntegral i * deltaTheta | i <- [0 .. numOrientation - 1]]
  --         weight
  --         pos
  --     idx = [(i, j) | i <- [-10 .. 10], j <- [-10 .. 10]]
  --     xs =
  --       L.map
  --         (\(i, j) ->
  --            let x = i * delta
  --                y = j * delta
  --             in ( (x, y)
  --                , L.sum $
  --                  cornerDistribution
  --                    delta
  --                    initScale
  --                    thetaSigma
  --                    tau
  --                    threshold
  --                    [ fromIntegral o * deltaTheta
  --                    | o <- [0 .. numOrientation - 1]
  --                    ]
  --                    weight
  --                    (x, y)))
  --         idx
  -- print . L.sortOn snd $ xs
  -- printf "%f\n" (L.sum cornerDist)
  -- print .
  --   L.sortOn snd .
  --   L.zip
  --     [fromIntegral i * deltaTheta / pi * 180 | i <- [0 .. numOrientation - 1]] $
  --   cornerDist
  -- plotImageRepa (folderPath </> "test.png") .
  --   ImageRepa 8 .
  --   fromListUnboxed (Z :. (1 :: Int) :. (21 :: Int) :. (21 :: Int)) .
  --   snd . L.unzip $
  --   xs
  -- let maxIdx =
  --       fst .
  --       L.maximumBy (\(_, a) (_, b) -> compare a b) .
  --       L.zip
  --         [ fromIntegral i * deltaTheta / pi * 180
  --         | i <- [0 .. numOrientation - 1]
  --         ] $
  --       cornerDist
  -- maxDist <- cornerDistribution' delta initScale thetaSigma tau threshold [fromIntegral i * deltaTheta | i <- [0 .. numOrientation - 1]] weight pos (maxIdx / 180 * pi)
  -- print maxIdx
  -- print maxDist
  flag <- doesFileExist histFilePath
  initialise []
  devs <- M.mapM device deviceIDs
  ctxs <- M.mapM (\dev -> CUDA.create dev []) devs
  ptxs <- M.mapM createTargetFromContext ctxs
  hist <-
    if flag
      then do
        printCurrentTime "read from files..."
        decodeFile histFilePath
      else sampleCartesianCorner
             histFilePath
             folderPath
             ptxs
             numPoints
             periodEnv
             delta
             numOrientation
             initScale
             thetaSigma
             tau
             threshold
             s
             phiFreq
             rhoFreq
             thetaFreq
             scaleFreq
             stdR2
             weight
  printCurrentTime "Done"
  printCurrentTime "Start Convloution.."
  plan <-
    makePlan
      folderPath
      emptyPlan
      numPointsRecon
      numPointsRecon
      numR2Freq
      (2 * thetaFreq + 1)
      (2 * scaleFreq + 1)
      (2 * phiFreq + 1)
      (2 * rhoFreq + 1)
  printCurrentTime "Compute DFT Plan done."
  let coefficients = getNormalizedHistogramArr hist
  harmonicsArray <-
    createHarmonics
      numR2Freq
      phiFreq
      rhoFreq
      thetaFreq
      scaleFreq
      (-s)
      periodR2
      periodEnv
      coefficients
  let initDistSource =
        computeInitialDistributionFourierPinwheel
          numR2Freq
          periodR2
          periodEnv
          phiFreq
          rhoFreq
          thetaFreq
          scaleFreq
          initSource
      initDistSink =
        computeInitialDistributionFourierPinwheel
          numR2Freq
          periodR2
          periodEnv
          phiFreq
          rhoFreq
          thetaFreq
          scaleFreq
          initSink
  source <- convolve harmonicsArray initDistSource
  sink <- convolve harmonicsArray initDistSink
  sourceR2 <- plotFPArray plan (folderPath </> "Source.png") source
  sinkR2 <- plotFPArray plan (folderPath </> "Sink.png") sink
  plotRThetaDist
    (folderPath </> "Source_Theta.png")
    (folderPath </> "Source_R.png")
    numPointsRecon
    360
    90
    periodEnv
    pos
    sourceR2
  plotRThetaDist
    (folderPath </> "Sink_Theta.png")
    (folderPath </> "Sink_R.png")
    numPointsRecon
    360
    90
    periodEnv
    pos
    sinkR2
  sinkR2TR <-
    computeUnboxedP .
    timeReversalRepa [fromIntegral (-thetaFreq) .. fromIntegral thetaFreq] $
    sinkR2
  completionR2 <- completionFieldRepa plan sourceR2 sinkR2TR
  (sumP . sumS . R.map (\x -> magnitude x ** 2) . rotate4D2 $ completionR2) >>=
    plotImageRepa (folderPath </> "Completion.png") .
    ImageRepa 8 .
    computeS .
    extend (Z :. (1 :: Int) :. All :. All) -- . R.map (\x -> log (x + 1))

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE QuantifiedConstraints #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
module Q.Interpolation where
import qualified Q.SortedVector as SV
import Numeric.GSL.Interpolation
import qualified Numeric.LinearAlgebra as V (Vector, fromList)
import Foreign (Storable)
import Data.List
class (Ord k, Storable k, Storable v) => Interpolator a k v where
  interpolate :: a -> [(k, v)] -> k -> v

class (Ord k, Storable k, Storable v) => InterpolatorV a k v where
  interpolateV :: a -> SV.SortedVector k -> V.Vector v -> k -> v

instance (Ord k, Storable k, Storable v, InterpolatorV a k v) => Interpolator a k v where
  interpolate a pts = interpolateV a xs' ys' where
    (xs, ys) = (unzip . sortOn fst) pts
    xs'      = SV.fromSortedList xs
    ys'      = V.fromList ys

max_stars_repo_name: ghais/lowq

max_issues_repo_name: ghais/lowq

max_forks_repo_name: ghais/lowq


text:
module Lib.Practice.P5
    (
        f5
    ) where

import Numeric.LinearAlgebra
import Lib.Utility

f5 = do
    let sigmoid x = 1 / (1 + exp(-x))
    let calcA x w b = add b $ x <> w
    let x = (1><2) [1,0.5 :: Double]
    let w1 = (2><3) [0.1,0.3,0.5,0.2,0.4,0.6]
    let b1 = (1><3) [0.1,0.2,0.3]
    let a1 = calcA x w1 b1
    let z1 = cmap sigmoid a1
    print z1

    let w2 = (3><2) [0.1,0.4,0.2,0.5,0.3,0.6]
    let b2 = (1><2) [0.1,0.2]
    let a2 = calcA z1 w2 b2
    let z2 = cmap sigmoid a2
    print z2

    let w3 = (2><2) [0.1,0.3,0.2,0.4]
    let b3 = (1><2) [0.1,0.2]
    let a3 = calcA z2 w3 b3
    let y = a3
    print y

    let a4 = (1><3) [1010,1000,990 :: Double]
    print $ softmax' a4

softmax' :: Matrix Double -> Matrix Double
softmax' a = expA / sumExpA
    where (r, c) = size a
          maxs = col $ map maxElement $ toRows a
          expA = cmap exp $ a - maxs
          sumExpA = row $ map sumElements $ toRows expA