text:
module Window (

    hamming
  , hammingC
  , c
  , Window
  , CWindow

  ) where

import Data.Vector
import qualified Data.Complex as C

type Window = Vector Double -- Windowing function

type CWindow = Vector (C.Complex Double) -- Windowing function

-- quick conversion to complex
c :: Num a => a -> C.Complex a
c = (C.:+ 0)

-- A hamming window of size n.
hamming :: Int -> Window
hamming m = generate m hamming' where
  hamming' n = 0.54 - 0.46*cos(2 * pi * (fromIntegral n)/(fromIntegral m-1))

-- A hamming window of type Complex Double of size n
hammingC :: Int -> CWindow
hammingC m = generate m hamming' where
  hamming' n = c $ 0.54 - 0.46*cos(2 * pi * (fromIntegral n)/(fromIntegral m-1))

ext: hs

max_stars_repo_name: davlum/haskell-stft

max_issues_repo_name: davlum/haskell-stft

max_forks_repo_name: davlum/haskell-stft


text:
module Main where

import Mandelbrot
import Data.Complex (Complex((:+)))
import Graphics.UI.SDL.Video
import Graphics.UI.SDL.Event (quitRequested)
import Graphics.UI.SDL.Types (Surface(surfacePixels))
import Foreign.C.String (newCString)
import Foreign.C.Types (CInt(..))
import Foreign.Ptr
import Foreign.Storable
import Data.Int (Int32(..), Int64(..))
import Data.Word
import Control.Monad
import Control.Concurrent (threadDelay)
import Data.Array
import Data.Array.MArray
import Data.Array.IO
import Data.IORef
import System.Random
import System.IO.Unsafe

sdlWindowPosCentered :: CInt
sdlWindowPosCentered = 0x2fff0000

iterations :: Int
iterations = 100

width :: Int
width = 1024

height :: Int
height = 786

widthF :: Double
widthF = fromIntegral width

heightF :: Double
heightF = fromIntegral height

mandelbrotI :: RealFloat a => Complex a -> Int
mandelbrotI = mandelbrot iterations
{-# INLINE mandelbrotI #-}

while :: Monad m => m Bool -> m a -> m ()
while test action = do
    val <- test
    if val
    then action >> while test action
    else return ()

main :: IO ()
main = do
    videoInit nullPtr
    title <- (newCString "Handelbrot")

    palette <- array (1, iterations) <$> sequence ( [sequence (i, (fromIntegral <$> randomRIO (0 :: Word, 0xffffff)) :: IO Word32) | i <- [1..iterations-1]] ++ [return (iterations, 0)] )

    window <- createWindow title sdlWindowPosCentered sdlWindowPosCentered (CInt $ fromIntegral width) (CInt $ fromIntegral height) 0
    surface <- getWindowSurface window
    lockSurface surface
    pixels <- liftM (castPtr.surfacePixels) $ peek surface :: IO (Ptr Word32)

    forM_ [0..height-1] (\y_ -> do
        forM_ [0..width-1] (\x_ -> do
            let x = 3/widthF * (fromIntegral x_) - 2
                y = -2/heightF * (fromIntegral y_) + 1
                mand = mandelbrotI (x :+ y)
            pokeElemOff pixels (x_ + y_*width) $ palette ! mand)
        updateWindowSurface window)

    unlockSurface surface
    updateWindowSurface window

    while (liftM not $ quitRequested) $ threadDelay 100000

max_stars_repo_name: hator/handelbrot

max_issues_repo_name: hator/handelbrot

max_forks_repo_name: hator/handelbrot


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GADTs     #-}
{-|
Module      : Grenade.Utils.PascalVoc
Description : Post-processing for TinyYoloV2
-}

module Grenade.Utils.PascalVoc (DetectedObject, labels, processOutput) where

import           Data.List                    (elemIndices)
import qualified Numeric.LinearAlgebra.Data   as NLA
import qualified Numeric.LinearAlgebra.Static as H

import           Grenade.Core.Shape
import           Grenade.Types

-- | List of possible labels of TinyYoloV2
labels :: [String]
labels
  = [ "aeroplane"
    , "bicycle"
    , "bird"
    , "boat"
    , "bottle"
    , "bus"
    , "car"
    , "cat"
    , "chair"
    , "cow"
    , "diningtable"
    , "dog"
    , "horse"
    , "motorbike"
    , "person"
    , "pottedplant"
    , "sheep"
    , "sofa"
    , "train"
    , "tvmonitor"
    ]

-- | Type alias for output from Yolo neural network.
--
--   Stores (x position, y position, width height, confidence score, object label).
type DetectedObject = (Int, Int, Int, Int, RealNum, String)

-- | Given the output of TinyYoloV2, finds all bounding boxes
--   with confidence higher than a threshold, and their labels
processOutput :: S ('D3 13 13 125) -- ^ Output tensor from TinyYoloV2
              -> RealNum           -- ^ Minimum confidence for bounding boxes, between 0 and 1. A good default is 0.3.
              -> RealNum           -- ^ Maximum IoU ratio to consider boxes not overlapping, between 0 and 1. A good default is 0.5.
              -> [DetectedObject]
processOutput (S3D mat) confThreshold iouThreshold = collapseBoxes iouThreshold $ map toDetectedObject filtered
  where
    filtered = filter (\(_, _, _, _, c, probs) -> snd (argMax probs) * c > confThreshold ) boxDescs
    out      = H.extract mat
    boxDescs = [ getBoundingBoxDesc out i j box | i <- [0..12], j <- [0..12], box <- [0..4]]

    toDetectedObject (x, y, w, h, c, probs) = (l, r, t, b, c, labels !! (fst $ argMax probs))
      where
        l = floor $ x - w / 2
        r = floor $ x + w / 2
        t = floor $ y - h / 2
        b = floor $ y + h / 2

    -- Returns (index of maximum element, maximum element)
    argMax probs = (i, m)
      where
        m = maximum probs
        i = head $ elemIndices m probs

-- | Given an IoU threshold and a list of boxes, find overlapping boxes and collapse them into one.
--
--   Boxes are determined to be overlapping if their intersection over union value passes the threshold.
--   When boxes are merged we pick the box with higher confidence score.
--
--   IoU threshold should be between 0 and 1, we've found 0.5 to work quite well.
collapseBoxes :: RealNum -> [DetectedObject] -> [DetectedObject]
collapseBoxes threshold xs = filter (\y -> all (doesNotOverlapWith y) xs) xs
  where
    -- Consider objects as overlapping if they have same label and their intersection over union area is too big.
    -- Only consider overlapping the object which has lower confidence.
    doesNotOverlapWith :: DetectedObject -> DetectedObject -> Bool
    doesNotOverlapWith y@(left, right, top, bottom, conf, label) x@(left', right', top', bottom', conf', label')
      = label /= label' || y == x || iou < threshold || conf > conf'
        where
          ix1   = max left left'
          ix2   = min right right'
          iy1   = max top top'
          iy2   = min bottom bottom'
          area  = (right  - left)  * (bottom  - top)
          area' = (right' - left') * (bottom' - top')
          inter = if ix1 < ix2 && iy1 < iy2 then (ix2 - ix1) * (iy2 - iy1) else 0
          iou   = fromIntegral inter / fromIntegral (area + area' - inter)

-- | Given the output matrix from YOLO, a cell position cy, cx and a bounding box b,
--   outputs (x, y, width, height, confidence score, probability distribution over 20 classes)
--   for said box.
--
--   Must have that 0 <= cy, cx <= 12 and 0 <= b <= 4.
getBoundingBoxDesc :: NLA.Matrix RealNum -> Int -> Int -> Int -> (RealNum, RealNum, RealNum, RealNum, RealNum, [RealNum])
getBoundingBoxDesc out cy cx b = (x, y, w, h, confidence * (maximum classes), classes)
  where
    anchors    = [1.08, 1.19, 3.42, 4.41, 6.63, 11.38, 9.42, 5.11, 16.62, 10.52]
    numClasses = 20
    channel    = b * 25

    tx = NLA.atIndex out ((channel + 0) * 13 + cy, cx)
    ty = NLA.atIndex out ((channel + 1) * 13 + cy, cx)
    tw = NLA.atIndex out ((channel + 2) * 13 + cy, cx)
    th = NLA.atIndex out ((channel + 3) * 13 + cy, cx)
    tc = NLA.atIndex out ((channel + 4) * 13 + cy, cx)

    x = (fromIntegral cx + sigmoid tx) * 32
    y = (fromIntegral cy + sigmoid ty) * 32

    w = exp tw * (anchors !! (2 * b    )) * 32
    h = exp th * (anchors !! (2 * b + 1)) * 32

    confidence = sigmoid tc

    classes = softmax [ NLA.atIndex out ((channel + 5 + c) * 13 + cy, cx) | c <- [0..numClasses - 1] ]

    sigmoid :: RealNum -> RealNum
    sigmoid = (1 /) . (1 +) . exp . negate

    softmax :: [RealNum] -> [RealNum]
    softmax xs = map (/sumxs) expxs
      where
        expxs = map exp xs
        sumxs = sum expxs

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
import           Advection1D
import           Basis
import           Element
import           Mesh
import           Node
import           Numeric.LinearAlgebra.HMatrix as HMat
import           Projection
import           ShapeFcns

linBasis   = Lagrange 1
linShpFcn  = TensorProduct linBasis 1
linShpFcn2 = TensorProduct linBasis 2
nodes      = [Node 0 [0.0], Node 1 [0.25]]
nodes2d    = [Node 0 [-1.0, -1.0], Node 1 [1.0, -1.0], Node 2 [1.0, 1.0], Node 3 [-1.0, 1.0]]
lineElem   = StructElem nodes 0
quadElem   = StructElem nodes2d 0
someCoords = [0.0, 0.0]
mesh2D     = generateMesh [0.0, 0.0] [1.0, 1.0] [3, 3] StructElem
mesh1D     = generateMesh [0.0] [1.0] [6] StructElem
meshDiscon = genDMesh (0.0::Double) (3.0::Double) 4 StructElem
lineElemMats = elemMatrices lineElem linShpFcn 1
globalMats = assembleGlobalMatrices meshDiscon linShpFcn 1 1.0

someFunction :: [Double] -> Double
someFunction [x] = sin x

nodalInitialConds = map (\e -> projectL2 e linShpFcn 2 someFunction) (getMeshElements meshDiscon) -- Should be [Maybe (Matrix Double)]
nodeCs = map (concatMap nodeCoordinates . getElementNodes) (getMeshElements meshDiscon)
dispElemIC :: [Double] -> Maybe (Matrix Double) -> String
dispElemIC nodeCoords (Just nodeICs) = concat [show (nodeCoords !! i) ++ " " ++ show (atIndex nodeICs (i,0)) ++ "\n" | i <- [0..1]]
dispElemIC _ Nothing = error "Linear solve failed."
icStrings = concat (zipWith dispElemIC nodeCs nodalInitialConds)

main = do
  print $ computeJacobian lineElem linShpFcn [0.0 :: Double]
  print $ show quadElem
  print $ show $ boundaryElements mesh2D !! 2
  print $ map (fst . elemMatrices lineElem linShpFcn) [1,2,3]
  writeMesh meshDiscon "mesh.dat"
  print $ HMat.toDense $ head globalMats
  print $ HMat.toDense $ globalMats !! 1
  print $ HMat.toDense $ last globalMats
  writeFile "sine_ic.dat" icStrings


text:
module Test.Test where

import Test.Framework (defaultMainWithArgs, testGroup, Test)
import Test.Framework.Providers.HUnit (testCase)
import Test.Framework.Providers.QuickCheck2 (testProperty)

import Test.QuickCheck (Property, NonNegative(..))
import Test.QuickCheck.Monadic (assert, monadicIO, run)
import Test.HUnit (assertBool)

import qualified Numeric.LinearAlgebra as LA
--import qualified Graphics.Rendering.Plot.HMatrix as G

import qualified Math.FROG.Tools as FR
import Math.FROG.Retrieval
import Math.FROG.Types
--import Test.Types

type PulseFn = Double -> Double -> LA.Complex Double

-----------------
-- Example Pulses
-----------------
singleGaus :: PulseFn
singleGaus sz k = exp ((-pi) / sz * (k - sz/2.0)^2) LA.:+ 0.0

doubleGaus :: PulseFn
doubleGaus sz k = (singleGaus sz (k - (sz / 8.0))) + (singleGaus sz (k + (sz / 8.0)))

mkPulse :: PulseFn -> Int -> ComplexSignal
mkPulse form sz = LA.buildVector sz mkElem
    where
    dsz = (fromIntegral sz) :: Double
    mkElem = form dsz . fromIntegral

caseRetrieval :: ComplexSignal -> IO ()
caseRetrieval field = do
    let sz = LA.dim field
    let dsz = (fromIntegral sz) :: Double
    let trc = FR.mkTrace field field

    o <- retrieve SHG trc

    let rfield = FR.center . FR.flatPolar2Complex $ o
    let mags = (LA.mapVector LA.magnitude) rfield

    putStrLn . show $ calcLoss SHG trc o

    --G.mplot [LA.linspace sz (0::Double, dsz-1), mags]

    --G.imshow trc

    --G.imshow (FR.mkTrace rfield rfield)

-----------------
---- List of Tests
-----------------

tests :: [Test]
tests =
    [ testGroup "Retrieval" [ testCase "Gaussian" $ caseRetrieval (mkPulse singleGaus 32)
                            , testCase "Double Gaussian" $ caseRetrieval (mkPulse doubleGaus 32)
                            ]
    ]

-- convenient whilst in ghci
runAllTests :: IO ()
runAllTests = defaultMainWithArgs tests []

runTests :: String -> IO ()
runTests p = defaultMainWithArgs tests ["--select-tests", p]

runGroup :: Int -> IO ()
runGroup i = defaultMainWithArgs [tests !! i] []

max_stars_repo_name: leroix/haskell-spsa-frog

max_issues_repo_name: leroix/haskell-spsa-frog

max_forks_repo_name: leroix/haskell-spsa-frog


text:
module Cube where

import Prelude hiding (Left, Right)
import Data.Function (on)
import Numeric.LinearAlgebra (Z, Vector)
import qualified Data.List as List
import qualified Numeric.LinearAlgebra as LA

type SideField = (Int, Int)

fieldsOnSide :: [SideField]
fieldsOnSide = [(x,y)|x <- [0..2], y <- [0..2]]

data Side = Top | Front | Right | Back | Left | Bottom
  deriving (Show, Read, Eq, Enum, Ord, Bounded)

data Color = White | Green | Red | Blue | Orange | Yellow
  deriving (Show, Read, Eq, Enum, Ord, Bounded)

type Cube = [(Side, Int, Int, Color)]

solvedCube :: Cube
solvedCube = do
  (side, color) <- zip [Top ..] [White ..]
  (x, y) <- fieldsOnSide
  [(side, x, y, color)]

computeIndex :: (Side, Int, Int, Color) -> Int
computeIndex (_, x, y, c) = x * 3 + y + fromEnum c * 9

cubeToVector :: Cube -> Vector Z
cubeToVector = LA.fromList . fmap go . List.sortBy (compare `on` project)
  where
    project (a,b,c,_) = (a,b,c)

    go :: (Side, Int, Int, Color) -> Z
    go (_, _, _, c) = fromIntegral $ fromEnum c

vectorToCube :: Vector Z -> Cube
vectorToCube = zipWith go solvedCube . fmap toColor . LA.toList
  where
    toColor :: Z -> Color
    toColor = toEnum . fromIntegral

    go :: (Side, Int, Int, Color) -> Color -> (Side, Int, Int, Color)
    go (s, x, y, _) c = (s, x, y, c)

oppositeSides :: Side -> Side -> Bool
oppositeSides Top Bottom = True
oppositeSides Front Back = True
oppositeSides Right Left = True
oppositeSides x y
  | fromEnum x > fromEnum y = oppositeSides y x
  | otherwise = False

adjacentSides :: Side -> Side -> Bool
adjacentSides x y
  | x == y = False
  | otherwise = not $ oppositeSides x y

max_stars_repo_name: runjak/hRubiks

max_issues_repo_name: runjak/hRubiks

max_forks_repo_name: runjak/hRubiks


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module JuliaC where

import Control.Monad (replicateM_)
import Control.Monad.Trans (lift)
import Control.Monad.Trans.State.Strict

import Data.Complex
import Data.Foldable (maximumBy)
import Data.Function (on)

import Firestarter

import Linear

import Numeric.AD

data Julia = Julia (V2 Double) Double

instance Metric Complex

instance Firestarter Julia V2 Double where
  isInside (Julia _ e) v = juliaDistance v <= e
  distanceAndDirection (Julia _ _) v =
    let
      (d,dir) = grad' juliaDistance v
      n = fmap signum $ maximumBy (compare `on` norm) $ map (negate dir*) basis
    in
      (d,n)
  startingLocation _ = V2 2 2

juliaFunction
  :: (Floating a, RealFloat a, Ord a)
  => Complex a
  -> StateT (Complex a,Complex a) (Either (Complex a, Complex a)) ()
juliaFunction c = do
  (q,dq) <- get
  if quadrance q > 200
    then lift $ Left (q,dq)
    else put (q*q+c, 2*^ (q*dq))

juliaDistanceEstimate :: Floating a => (Complex a, Complex a) -> a
juliaDistanceEstimate (q,dq) = 0.5 * r * log r / dr
  where r = norm q
        dr = norm dq

juliaDistance :: (Floating a, Ord a, RealFloat a) => V2 a -> a
juliaDistance v = juliaDistanceEstimate . either id id $ execStateT s qs
  where
    s = replicateM_ 100 $ juliaFunction c
    c = initC
    qs = fromVec v

initC :: Floating a => Complex a
initC = (-0.125) :+ (-0.256)
--initComplex = (-0.4) :+ 0.6

fromVec :: Floating a => V2 a -> (Complex a, Complex a)
fromVec (V2 x y) = (x :+ y, 1 :+ 0)

max_stars_repo_name: fabricatedmath/firestarter

max_issues_repo_name: fabricatedmath/firestarter

max_forks_repo_name: fabricatedmath/firestarter


text:
{-# OPTIONS_HADDOCK not-home #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE ViewPatterns #-}
module Types
    ( Env, GlobalEnv, LocalEnv
    , Arity (..)
    , InTail
    , Builtin
    , Primitive (..)
    , Macro (..)
    , Val (..), Symbol, Port(..)
    , RealNumber (..), Number (..)
    , symbolName, symbolAsString
    , Ref
    , LispErr (..), ExceptionCont (..)
    , ThrowsError
    , IOThrowsError
    , extractValue
    , isTerminationError
    , liftEither
   -- , runIOThrows
    , EM (..), emCont
    , EvalState (..), dynPoint, currentHandlers
    , StackFrame (..)
    , DynamicEnv(DynamicEnv)
    , ExceptionHandler
    , DynamicPoint (..)
    , Opts
   -- , TraceType
    , liftIOThrows
    ) where

import Control.Monad.Cont
    ( (<=<),
      MonadIO(..),
      MonadCont(..),
      ContT(ContT),
      cont,
      runCont,
      Cont )
import Control.Monad.Except           ( ExceptT, MonadError(..), liftEither
                                      , runExceptT )
import Control.Monad.State.Lazy       ( MonadState(get, put), gets, modify )
import Data.Complex                   ( Complex(..) )
import Data.HashMap.Strict            ( HashMap )
import Data.IORef                     ( IORef, readIORef, writeIORef )
import Data.Ratio                     ( (%), denominator, numerator )
import Data.Text                      ( Text, unpack )
import qualified Data.Text           as T
import qualified Data.Vector         as V
import qualified Data.Vector.Mutable as V
import qualified Data.Vector.Unboxed as U
import           Data.Word                      ( Word8 )
import           System.IO                      ( Handle )

import Options

{- Note: [IORefs in Envs]

Vals already have a notion of being "primitive" or a "pointer".
So why do 'Env's store 'IORef's? We want to be able to
_duplicate environments_ (for the purposes of capturing them in closures)
and be sure that the copied environments reflect changes in the originals.

It's not impossible that I'm completely off-base here and duplicating the
environments is not actually worth it in terms of performance. A duplicated
environment is always smaller than the original, since it flattens the
scope structure and removes shadowed bindings. However, That _usually_
isn't significantly smaller, and currently we pay an extra dereference for
every single name lookup. Low-hanging fruit?
-}
type Env = IORef (HashMap Symbol (IORef Val))
type GlobalEnv = Env
type LocalEnv  = [Env] -- ^ stack of environments forms local env
type Ref = IORef

-- * Function types and components
--type Arity = Int
data Arity
  = Exactly Int
  | AtLeast Int
  | Between Int Int
  deriving (Eq, Ord)
type InTail = Bool
type Builtin = [Val] -> EM Val
data Primitive = Prim Symbol Arity Builtin
-- The additional 'InTail' flag is necessary for macros like 'begin'
-- and 'if' which need to use a different 'eval' behavior if they are
-- in tail position. We do global CPS anyway (EM monad) but this information
-- is needed to maintain the debug stack.
data Macro = Macro Arity (InTail -> Builtin)

-- | Scheme values.
--
-- Some Scheme values can be either mutable or immutable. These are strings,
-- lists, and vectors. The immutable variant begins with I, where applicable.
--
-- The mutable data structures are interpreter-only;
-- the parser cannot construct them.
--
-- Invariant: All objects contained in an IPair or IVector are immutable.
--   Note that this is implicitly transitive.
data Val
  = Symbol Symbol

  | Number !Number

  | String !(Ref Text)
  | IString Text
  | Char Char
  | Bool Bool

  | Primitive Arity Builtin Symbol
  | PrimMacro Arity (InTail -> Builtin) Symbol
  | Continuation DynamicPoint (Val -> EM Val)
  | Closure
     { params :: [Symbol]
     , vararg :: Maybe Symbol
     , body   :: [Val]
     , cloEnv :: LocalEnv
     , name   :: Maybe Symbol
     }

  | Port Port

  | Nil
  | Pair !(Ref Val) !(Ref Val)
  -- this can't be made strict, else, attempting to construct an
  -- cyclic static data structure containing a pair /anywhere in the structure/
  -- would hang!
  | IPair Val Val

  | Vector  !(V.IOVector Val)
  | IVector !(V.Vector Val)
  | IByteVector !(U.Vector Word8)

   -- if we ever add a way to register writers for custom
   -- record types, then error objects should become a
   -- standard library record type and this should go.
  | Error Text [Val]
  | Exception LispErr
  | EOF

    -- | Used to identify when a recursive binding refers to itself during
    -- initialization, which is an error.
    -- If 'unspecified' or 'void' are added in the future, they should be
    -- distinct from this.
  | Undefined
    -- | Implements multiple values for the purposes of call-with-values.
    -- If one of these ever appears somewhere that it shouldn't, it will
    -- eventually get sent through `eval`, which will raise an error.
  | MultipleValues [Val]

type Symbol = Text

symbolName :: Symbol -> Text
symbolName (T.uncons -> Just ('~', name)) = name
symbolName (T.uncons -> Just ('&', name)) =
  let pieces = T.splitOn "~" name
      -- pieces are some number of components from the name
      -- if it contained tildes, followed by two components
      -- which are the unique identifiers.
  in T.intercalate "~" $ reverse $ drop 2 $ reverse pieces
symbolName primName = primName

-- | Converts a symbol to its human-readable name as a string.
-- Mostly useful for Show/write and in 'panic' messages.
symbolAsString :: Symbol -> String
symbolAsString = unpack . symbolName

data Port
  = HandlePort !Handle
  -- consuming is permanent, as indexing characters
  -- into a Unicode stream is necessarily O(n)
  -- so supporting seeking is horribly slow.
  -- When closed, the IORef is cleared so that the
  -- text can eventually be freed.
  | InputTextPort !(IORef Text) !(IORef Bool) -- open?
  -- collect chunks, which will later be collected
  -- by 'T.concat'. The chunks are stored in LIFO order.
  -- When closed, the IORef is cleared so that the
  -- list can be reclaimed by GC.
  | OutputTextPort !(IORef [Text]) !(IORef Bool) -- open?
  deriving (Eq)

-- | Scheme numbers. Note that there is no Ord Number instance;
-- complex numbers are not orderable and Scheme comparisons of
-- complex numbers are errors.
data Number
  = Real {-# UNPACK #-}!RealNumber
  | Complex {-# UNPACK #-}!(Complex RealNumber)

data RealNumber
  = Bignum {-# UNPACK #-}!Integer
  | Flonum {-# UNPACK #-}!Double
  | Ratnum {-# UNPACK #-}!Rational

-- | This is equality in the sense of Scheme '='.
instance Eq RealNumber where
  Flonum f == Flonum g = f == g
  Flonum f == bigOrRat = not (isInfinite f) && r1 == r2
    -- note that we make the fp number exact instead of making the
    -- exact number fp. This is critical. An example from the R7RS
    -- paper is: let big be (exp 2 1000). Assume big is exact and
    -- that inexact numbers are IEEE doubles [they are]. Then
    -- (= (inexact big) (+ big 1)) should be #f, but if inexact
    -- conversion is used, it will be #t.
    where r1 = toRational f
          r2 = case bigOrRat of
            Bignum i -> i % 1
            Ratnum r -> r
            Flonum{} -> error "(==)@RealNumber: impossible flonum"
  bigOrRat == n@Flonum{} = n == bigOrRat
  Bignum i == Bignum j = i == j
  Bignum i == Ratnum r = i % 1 == r
  Ratnum r == Bignum i = r == i % 1
  Ratnum r == Ratnum s = r == s

instance Eq Number where
  Real a == Real b = a == b
  Real a == Complex (r :+ i) = a == r && i == Bignum 0
  Complex (r :+ i) == Real b = r == b && i == Bignum 0
  Complex (a :+ b) == Complex (x :+ y) = a == x && b == y

instance Show RealNumber where
  show (Bignum i) = show i
  show (Flonum f)
    | isInfinite f = sign : "inf.0"
    | isNaN f      = "+nan.0"
    | otherwise    = show f
    where sign = if f < 0 then '-' else '+'
  show (Ratnum r)
    | denominator r == 1 = show (numerator r)
    | otherwise = show (numerator r) ++ "/" ++ show (denominator r)

-- | Doesn't compare mutable things; most likely you want
--
-- (1) 'eqSSH':    referential equality
-- (2) 'eqvSSH':   value equality for atoms, otherwise referential
-- (3) 'equalSSH': structural equality
eqVal :: Val -> Val -> Bool
eqVal (Symbol s) (Symbol s') = s == s'
eqVal (Number n) (Number n') = n == n'
eqVal (IString s) (IString s') = s == s'
eqVal (Char c) (Char c') = c == c'
eqVal (Bool b) (Bool b') = b == b'
eqVal (IVector v) (IVector v') = v == v'
eqVal (Port p) (Port p') = p == p'
eqVal (Primitive _ _ n) (Primitive _ _ n') = n == n'
eqVal Continuation{} _ = False -- perhaps it's possible to give continuations
eqVal _ Continuation{} = False -- unique identifiers to make this work.
                               -- However... why? 'eq?' works fine by just
                               -- comparing their StableNames and anything we
                               -- might want to test can be tested with 'eq?'.
eqVal (Closure p v b _ n) (Closure p' v' b' _ n') =
    and [ p == p'
        , v == v'
        , b == b'
        , n == n'
        ]
eqVal Nil Nil = True
eqVal p@IPair{} q@IPair{} =
    fromList p == fromList q
  where
    fromList :: Val -> ([Val], Val)
    fromList (IPair c d) =
      first (c:) $ fromList d
    fromList obj = ([], obj)

    first f (a, b) = (f a, b)
eqVal _ _ = False

-- Eq Val is defined here so that 'LispErr' can derive 'Eq'.
-- But it might be more sane to just not allow pure comparisons of
-- 'Val's or 'LispErr's. I think (?) it's only useful in tests
-- (and I'm not even sure that the tests use it anymore).

-- | Compares closures by shape rather than only by name/location tag.
-- Can't compare continuations (because they currently aren't tagged).
-- Can't compare mutable pairs (use 'valsSameShape' instead).
instance Eq Val where (==) = eqVal

-- | Only so that LispErr can derive Eq.
newtype ExceptionCont = EC (Val -> EM Val)
data LispErr = NumArgs Arity [Val]
             | TypeMismatch String Val
             | Parser String
             | BadSpecialForm Val
             | NotFunction String Val
             | UnboundVar String Symbol
             | EvaluateDuringInit Symbol
             | SetImmutable String -- type name
             | CircularList
             | EmptyBody

             | Condition (Maybe ExceptionCont) -- ^ continuable?
                         Val

             | Default String
             | Quit

-- This is not derived _purely_ so that we can compare
-- symbols by symbolName in two cases. Annoying.
instance Eq LispErr where
  NumArgs a1 vs1        == NumArgs a2 vs2        = a1 == a2 && vs1 == vs2
  TypeMismatch s1 v1    == TypeMismatch s2 v2    = s1 == s2 && v1 == v2
  Parser s1             == Parser s2             = s1 == s2
  BadSpecialForm v1     == BadSpecialForm v2     = v1 == v2
  NotFunction s1 v1     == NotFunction s2 v2     = s1 == s2 && v1 == v2
  UnboundVar s1 n1      == UnboundVar s2 n2      = s1 == s2 && symbolName n1 == symbolName n2
  EvaluateDuringInit n1 == EvaluateDuringInit n2 = symbolName n1 == symbolName n2
  SetImmutable s1       == SetImmutable s2       = s1 == s2
  CircularList          == CircularList          = True
  EmptyBody             == EmptyBody             = True
  Condition _ v1        == Condition _ v2        = v1 == v2
  Default s1            == Default s2            = s1 == s2
  Quit == Quit = True
  _ == _ = False

instance Eq ExceptionCont where
  EC _ == EC _ = True -- let the comparison be by everything else

type ThrowsError = Either LispErr

--trapError :: ThrowsError String -> ThrowsError String
--trapError action = catchError action (return . show)

extractValue :: ThrowsError a -> a
extractValue (Right val) = val
extractValue Left{} = error "extractValue: Left"

isTerminationError :: LispErr -> Bool
isTerminationError Quit = True
isTerminationError _    = False

type IOThrowsError = ExceptT LispErr IO

--runIOThrows :: IOThrowsError String -> IO String
--runIOThrows action = extractValue . trapError <$> runExceptT action

{- Note: [EM return types]
Ideally, we'd like EM to be able to return any type, instead of just 'Val'.
MTL's Cont type handles this with a type variable; 'Cont r a'.

I tried lots of variations on EM that didn't require putting an 'r' type
everywhere. I couldn't make it work. The most promising attempt was
type EMCont r = EvalState -> IO (Either LispErr r)
newtype EMR r a = EM { unEM :: Cont (EMCont r) a }
  deriving...
type EM a = forall r. EMR r a

This almost works. In fact, if we didn't have continuations, it _would_
work. Alas, the need to be able to reify an EM computation and store it
into a Val required either making Val quantify over r (non-universally)
or making Continuation existential. The first doesn't solve the problem,
it just pushes the problem down onto Val. The second doesn't work because
we can't apply an existentially quantified continuation in 'apply', because
the 'r' in the type of 'apply' is universally quantified.

So instead, we force EM actions to produce a Val. This is OK, because we
can use a hack with IORefs to get the actual value we care about instead,
if needed.
-}

{- Note: [EM structure]
EM uses a little trick about continuation monads to embed state and errors,
instead of stacking StateT and ErrorT. The reason for this is that
ConT r (ErrorT e m) a doesn't (and in fact, can't!) provide a MonadError
instance. However, this type can.

Cont r a ~ (a -> r) -> r, so
Cont EMCont a ~ (a -> EvalState -> IO (...)) -> EvalState -> IO (...)
-}
type EMCont = EvalState -> IO (Either LispErr Val, EvalState)
-- | The Evaluation Monad
newtype EM a = EM { unEM :: Cont EMCont a }
  deriving (Functor, Applicative, Monad)

emCont :: ((a -> EMCont) -> EMCont) -> EM a
emCont = EM . cont

emThrow :: LispErr -> EM a
emThrow e = do
  handlersRef <- gets currentHandlers
  handlers    <- liftIO $ readIORef handlersRef
  handler     <- case handlers of
    [] -> error "emThrow: no handlers"
    [handler] -> pure handler -- leave initial handler alone
    (handler:cdrHandlers) ->
      liftIO $ writeIORef handlersRef cdrHandlers >> return handler
  _ <- handler e
  -- If the handler returns, we are to raise an exception with an unspecified
  -- relationship to e. For now, we simply re-raise e itself.
  -- This is moderate future-proofing. If the 'Error' objects ever became
  -- a primitive record type instead of a constructor of 'Val', we'd be in
  -- trouble here.
  emThrow e

-- | State is recovered from the point that 'emCatch' was invoked.
-- This function is for the MonadError EM instance and should not
-- be used in practice, as it doesn't know anything about
-- exception handlers.
emCatch :: EM a -> (LispErr -> EM a) -> EM a
emCatch (EM m) restore = emCont $ \k s -> do
    -- this leaves the catch frame around until the entire continuation
    -- chain, including the continuation of the catch itself, is done.
    -- That seems wrong, so perhaps this should be:
    {- mr <- runEM m (\x s -> pure (Right x, s)) a s
       case mr of
         (Left e, _s0) -> runEM (restore e) k a s
         (Right v, s') -> runEM (return v)  k a s'
    -}
    -- Edit, that doesn't work for type reasons but can work with anIORef.
    -- But why bother? We don't use this function.
    mr <- runCont m k s
    case mr of
        (Left e, _s0)  -> runCont (unEM (restore e)) k s
        (Right _v, _s) -> pure mr

emGet :: EM EvalState
emGet = emCont $ \k s -> k s s

emPut :: EvalState -> EM ()
emPut s = emCont $ \k _s -> k () s

emLiftIO :: IO a -> EM a
emLiftIO io = emCont $ \k s -> io >>= flip k s

-- | Invoking a continuation restores the state to when it was captured.
instance MonadCont EM where
  callCC f = emCont $ \k s ->
    runCont (unEM $ f (\x -> emCont $ \ _ _ -> k x s)) k s
  --{-# INLINE callCC #-}

instance MonadState EvalState EM where
  get = emGet
  --{-# INLINE get #-}
  put = emPut
  --{-# INLINE put #-}

-- | This instance made sense historically, but do the inability of the type
-- of 'catchError' to be instantiated by a function which uses 'Val'-specific
-- exception handlers, it no longer makes sense.
--
-- Consider it deprecated. However, when it is removed, 'EvaluationMonad' will
-- export a function named 'throwError' which behaves as a Haskell-level
-- version of raise.
instance MonadError LispErr EM where
  throwError = emThrow
  catchError = emCatch

instance MonadIO EM where
  liftIO = emLiftIO
  --{-# INLINE liftIO #-}

instance MonadFail EM where
  fail s = throwError . Default $
    "The following error occurred, please report a bug.\n" ++ s
  --{-# INLINE fail #-}

instance HasOpts EM where
  getOpts = gets options
  setOpts opts = modify $ \s -> s { options = opts }

liftIOThrows :: IOThrowsError a -> EM a
liftIOThrows = liftEither <=< liftIO . runExceptT

-- | The current state of evaluation
data EvalState = ES { globalEnv  :: GlobalEnv
                    , localEnv   :: LocalEnv
                    , stack      :: [StackFrame]
                    , dynEnv     :: DynamicEnv
                    , options    :: Opts
                    }

type ExceptionHandler = LispErr -> EM Val

-- | Haskell-level model of the dynamic Scheme environment.
-- Be aware that this is not a dynamic-extent? object
-- ala SRFI-154. It is the single representation of whichever
-- dynamic extent happens to be current at this point in Thcheme's
-- execution. It mutates as Thcheme executes.
data DynamicEnv = DynamicEnv
  { deDynPoint :: Ref DynamicPoint
  , deHandlers :: Ref [ExceptionHandler]
  }

dynPoint :: EvalState -> Ref DynamicPoint
dynPoint = deDynPoint . dynEnv

currentHandlers :: EvalState -> Ref [ExceptionHandler]
currentHandlers = deHandlers . dynEnv

-- | Associate a call in the callTrace with its captured local environment.
data StackFrame = StackFrame Val (Maybe LocalEnv)

-- | Dynamic points; this is the core of the implementation of dynamic-wind.
-- The details of the implementation are mainly sourced from
-- https://www.cs.hmc.edu/~fleck/envision/scheme48/meeting/node7.html
--
-- Each dynamic point stores the "before" and "after" thunks from the
-- dynamic-wind call, as well as the "parent" node, which is either the
-- root, or the dynamic point of the next innermost dynamic-wind call.
-- The root is the Point whose parent is the Sentinel, and its
-- before/after ref contains unspecified values that should not be called.
--
-- The implementation of dynamic-wind spans a couple of modules.
-- Here, continuation objects are defined as holding a DynamicPoint as
-- well as a function; this is the dynamic point of the dynamic extent in
-- which that continuation was captured. When a continuation is captured
-- by callWithCurrentContinuation in Misc.hs, it captures the current value of
-- dereferencing dynPoint as well. When a continuation is applied by
-- 'call' in Evaluation.hs, dynPoint is rerooted to the captured DynamicPoint
-- before the continuation's function itself is called. The implementation of
-- dynamicWind in Misc.hs contains two more rerootDynPoint calls that
-- cause the normal enter-and-exit calls to before and after.
data DynamicPoint = Sentinel | Point (Ref (Val, Val)) (Ref DynamicPoint)
  deriving Eq

max_stars_repo_name: JKTKops/Thcheme

max_issues_repo_name: JKTKops/Thcheme

max_forks_repo_name: JKTKops/Thcheme


text:
-- |
-- Module      : AOC.Challenge.Day11
-- Copyright   : (c) Justin Le 2018
-- License     : BSD3
--
-- Maintainer  : justin@jle.im
-- Stability   : experimental
-- Portability : non-portable
--
-- Day 11.  See "AOC.Solver" for the types used in this module!

module AOC.Challenge.Day11 (
    day11a
  , day11b
  ) where

import           AOC.Common                       (Point, meanVar)
import           AOC.Solver                       ((:~>)(..))
import           Control.DeepSeq                  (force)
import           Data.Foldable                    (maximumBy)
import           Data.Ix                          (range)
import           Data.Map                         (Map)
import           Data.Maybe                       (catMaybes)
import           Data.Ord                         (comparing)
import           Data.Profunctor                  (dimap)
import           Linear                           (V2(..))
import           Text.Read                        (readMaybe)
import qualified Control.Foldl                    as F
import qualified Data.Map                         as M
import qualified Data.Set                         as S
import qualified Statistics.Distribution          as D
import qualified Statistics.Distribution.Binomial as D
import qualified Statistics.Distribution.Normal   as D

powerLevel :: Int -> Point -> Int
powerLevel sid (V2 x y) = hun ((rid * y + sid) * rid) - 5
  where
    hun = (`mod` 10) . (`div` 100)
    rid = x + 10

findMaxThree :: Map Point Int -> Point
findMaxThree mp = fst
                . maximumBy (comparing snd)
                . map (\x -> (x, go x))
                $ range (V2 1 1, V2 298 298)
  where
    go p = sum [ mp M.! (p + shift)
               | shift <- range (V2 0 0, V2 2 2)
               ]

mkMap :: Int -> Map Point Int
mkMap i = M.fromSet (powerLevel i) . S.fromList $ range (V2 1 1, V2 300 300)

day11a :: Int :~> Point
day11a = MkSol
    { sParse = readMaybe
    , sShow  = \(V2 x y) -> show x ++ "," ++ show y
    , sSolve = Just . findMaxThree . mkMap
    }

day11b :: Int :~> (Point, Int)
day11b = MkSol
    { sParse = readMaybe
    , sShow  = \(V2 x y, s) -> show x ++ "," ++ show y ++ "," ++ show s
    , sSolve = Just . findMaxAny . mkMap
    }

findMaxAny :: Map Point Int -> (Point, Int)
findMaxAny mp = fst $ go 1
  where
    go n
       | goOn > 0.001 = maximumBy (comparing snd) [((pMax, n), oMax), go (n + 1)]
       | otherwise    = ((pMax, n), oMax)
                      -- & traceShow (n, oMax, goOn)
      where
        (pMax, oMax) = maximumBy (comparing snd)
            [ (p, fromIntegral (fromSAT sat p n))
            | p <- range (V2 1 1, V2 (300 - n + 1) (300 - n + 1))
            ]
        goOn = sum
            [ probGreaterThan oMax n'
            | n' <- [n + 1 .. 300]
            ]
    !sat = summedAreaTable mp
    σ :: Double
    σ = sqrt $ ((4 + 5)**2)/12   -- stdev of uniform distribution between -5 and 4
    probGreaterThan o n
        | prob2 == 0 = 0
        | otherwise  = probIn2
      where
        n' = fromIntegral n
        distr2  = D.normalDistr (-(n' ** 2) * 0.5) (n' * σ)
        prob2   = D.complCumulative distr2 o
        numIn2  = ( (300 `div` n) ^ (2 :: Int)    -- we compensate for dependence
                  + (300 - n + 1)^(2 :: Int)
                  ) `div` 2
        probIn2 = 1 - D.probability (D.binomial numIn2 prob2) 0
                    -- this is technically not a bernoulli process,
                    -- because our items are dependent. but we fudge this
                    -- by tweaking the number of trials

fromSAT :: Map Point Int -> Point -> Int -> Int
fromSAT sat (subtract (V2 1 1)->p) n = sum . catMaybes $
    [            M.lookup p            sat
    ,            M.lookup (p + V2 n n) sat
    , negate <$> M.lookup (p + V2 0 n) sat
    , negate <$> M.lookup (p + V2 n 0) sat
    ]

summedAreaTable :: Map Point Int -> Map Point Int
summedAreaTable mp = force sat
  where
    sat = M.mapWithKey go mp
    go p0 v = (+ v) . sum . catMaybes $
      [ negate <$> M.lookup (p0 - V2 1 1) sat
      ,            M.lookup (p0 - V2 1 0) sat
      ,            M.lookup (p0 - V2 0 1) sat
      ]

-- | Debug: find the variance of the map at every square size
_chunkyVars :: Map Point Int -> Map Int Double
_chunkyVars mp = flip M.fromSet (S.fromList [1..300]) $ \n ->
    F.fold (dimap fromIntegral snd meanVar)
      [ fromSAT sat p n
      | p <- range (V2 1 1, V2 (300 - n + 1) (300 - n + 1))
      ]
  where
    !sat = summedAreaTable mp


text:
{-# LANGUAGE OverloadedLists
           , TypeFamilies
           , FlexibleContexts #-}

module Utils.Quaternions where

import           Number.Quaternion hiding (normalize,norm)
import qualified Number.Quaternion as Q
import           Numeric.LinearAlgebra (Matrix, Vector, tr, cmap,(<>), inv,single,double)
import qualified Data.Vector.Storable as VS
import qualified Numeric.LinearAlgebra.Data as LD

import Utils (mat3ToMat4)

qvMul :: VS.Vector Float -> T Float -> T Float
qvMul v = quatConcat (0 +:: (VS.unsafeIndex v 0,VS.unsafeIndex v 1,VS.unsafeIndex v 2))

quatToLInv inertiaDiag quat =
  let orientation = quatToMat3 quat
  in tr orientation
   <> (single . inv . double) inertiaDiag
   <> orientation

rotateWithQuat :: T Float -> VS.Vector Float -> VS.Vector Float
rotateWithQuat quat [a,b,c] =
  (\(a,b,c) -> [a,b,c]) .
  imag .
  flip quatConcat (qInverse quat) $
  quatConcat quat (0 +:: (a,b,c))

addQuat q1 q2 =
  let r1 = real q1
      r2 = real q2
      i1 = imag q1
      i2 = imag q2
  in (r1 + r2) +:: sum3 i1 i2

quatToMatV :: T Float -> VS.Vector Float
quatToMatV q =
  let r = real q
      (i,j,k) = imag q
      r2 = r*r
      i2 = i*i;   j2 = j*j;   k2 = k*k;
      ri = 2*r*i; rj = 2*r*j; rk = 2*r*k;
      ij = 2*i*j; ki = 2*k*i; kj = 2*k*j;
  in [ r2+i2-j2-k2,       ij-rk,       ki+rj, 0
     ,       ij+rk, r2-i2+j2-k2,       kj-ri, 0
     ,       ki-rj,       kj+ri, r2-i2-j2+k2, 0
     ,           0,           0,           0, 1]

quatToMat = tr . LD.reshape 4 . quatToMatV

quatToMat3 :: T Float -> Matrix Float
quatToMat3 = LD.fromArray2D . toRotationMatrix

sum3 (a,b,c) (d,e,f) = (a+d,b+e,c+f)

scale3 n (a,b,c) = (n*a,n*b,n*c)

quatConcat q0 q1 =
  let w0 = real q0
      w1 = real q1
      v0 = imag q0
      v1 = imag q1
      wr = w1 * w0 - scalarProduct v1 v0
      vr = scale3 w1 v0 `sum3` scale3 w0 v1 `sum3` crossProduct v1 v0
  in wr +:: vr

qInverse :: T Float -> T Float
qInverse q =
  let w = real q
      v = (\(a,b,c) -> (negate a, negate b, negate c)) $ imag q
  in w +:: v

qZero :: T Float
qZero = 1 +:: (0,0,0)

max_stars_repo_name: Antystenes/CPG

max_issues_repo_name: Antystenes/CPG

max_forks_repo_name: Antystenes/CPG


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

module Data.Vector.Image.Matrix (
  toMatrix
, fromMatrix
, toMatrixList
, fromMatrixList
, toMatrixD
, fromMatrixD
) where

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

import Data.Vector.Image.Color.RGB
import Data.Vector.Image.Image

toMatrix :: (Storable a)
        => Image a
        -> Matrix a
toMatrix (Image w _h v) = reshape w v

fromMatrix :: (Element a)
        => Matrix a
        -> Image a
fromMatrix mat = Image (cols mat) (rows mat) (flatten mat)

toMatrixList :: (ColorSpace c a,Storable (c a),Storable a,Element a)
              => Image (c a)
              -> [Matrix a]
toMatrixList img = P.map toMatrix $ splitColor img

fromMatrixList :: (ColorSpace c a,Storable (c a),Storable a,Element a)
              => [Matrix a]
              -> Image (c a)
fromMatrixList mats = mergeColor $ P.map fromMatrix mats

toMatrixD :: (ColorSpace c a,Storable (c a),Storable a,Element a,Integral a)
         => Image (c a)
         -> [Matrix Double]
toMatrixD img = P.map toMatrix $ P.map toDouble $ splitColor img

fromMatrixD :: (ColorSpace c a,Storable (c a),Storable a,Element a,Integral a)
           => [Matrix Double]
           -> Image (c a)
fromMatrixD mats = mergeColor $ P.map fromDouble $ P.map fromMatrix mats

max_stars_repo_name: junjihashimoto/img-vector

max_issues_repo_name: junjihashimoto/img-vector

max_forks_repo_name: junjihashimoto/img-vector


text:
module Main (main) where

-- package の import
-- CSVの読み込み
import qualified    Text.CSV                                as CSV
import qualified    Data.Text                               as T
import              Data.Text
import              Data.Attoparsec.Text                    as DAT

-- 可視化
import              Graphics.Rendering.Chart.Easy            hiding ( (:<))
import              Graphics.Rendering.Chart.Backend.Cairo
import              Graphics.Rendering.Chart.Axis
import              Graphics.Rendering.Chart.Axis.Int
import              Graphics.Rendering.Chart.Grid

-- 重回帰
import Statistics.Regression

-- その他
import qualified    Control.Monad                           as CM
import qualified    Data.List                               as L
import qualified    Data.Vector.Unboxed                     as VU

-- | Doubleのパーサ
{-# INLINE parseDouble #-}
parseDouble :: Text -> Double
parseDouble tx = case DAT.parseOnly DAT.double tx of
        Right r   -> r
        Left  l   -> error $ "Error on parseDouble : " ++ show l

type Title      = String
type FileName   = String
type Label      = String

-- | 折れ線グラフの作成
plotLine :: (PlotValue a, PlotValue b) => [(Label, [(a, b)])] -> FileName -> Title -> IO ()
plotLine xs file title
    = toFile def (file ++ ".png") $ do
        layout_title .= title
        CM.forM_ xs $ \(title,ys) -> plot $ line title [ys]

plotBar :: [(String,[Double])] -> FileName -> Title -> [Label] -> IO ()
plotBar xs file title labels
    = toFile def (file ++ ".png") $ do
    layout_title .= title
    layout_title_style . font_size .= 10
    layout_x_axis . laxis_generate .= autoIndexAxis (L.map fst xs)
    plot $ fmap plotBars $ bars labels (addIndexes (L.map snd xs))

-- | 散布図の作成
plotGridPoints :: (PlotValue a, PlotValue b)
             => [[[(Label, [(a, b)])]]] -> FileName -> Title -> IO ()
plotGridPoints xs file name
    = let  ys   = aboveN $ (flip L.map) xs
                $ \ cols  -> besideN
                $ (flip L.map) cols -- columns
                $ \ series -> layoutToGrid
                $ execEC
                $ CM.forM_ series
                $ \(t,xs) -> plot $ points t xs
    in CM.void  $ renderableToFile def (file ++ ".png")
                $ fillBackground def
                $ gridToRenderable
                $ title `wideAbove` ys
  where
    title = setPickFn nullPickFn $ label ls HTA_Centre VTA_Centre name
    ls = def { _font_size   = 15 , _font_weight = FontWeightBold }

main = do
    -- CSVファイルの読み込み
    -- Dataフォルダを作成し,そこにデータを入れておきましょう
    xs <- CSV.parseCSVFromFile "Data/multiple_regression_test.csv" >>= \res
       -> case res of
            Right xs -> return $ L.map (L.map T.pack) xs
            Left  xs -> error $ "error on parseCSVFromFile: " ++ show xs

    -- Headerの抽出 Vectorへの変換
    let headers = L.map T.unpack
                $ L.head xs

    -- 値のDoubleへの変換 Vectorへの変換
    let body    = L.transpose
                $ L.map (L.map parseDouble)
                $ L.tail xs

    let year        = body !! 0
        death       = body !! 1
        index       = body !! 2
        hospital    = body !! 3

    -- データをグラフで表示
    plotLine    (L.map (\(x,y) -> (x, L.zip year y))
                [("Number of Death", death)
                ,("Economic Indicators", index)
                ,("Number of Hospitals", hospital)])
                "Data/initial_data"
                "Data Visualization"

    -- 散布図行列を作成する
    plotGridPoints  [[[(xt ++ "_" ++ yt, L.zip xs ys)]
                     | (yt,ys) <- L.zip (L.tail headers) (L.tail body)]
                    | (xt,xs) <- L.zip (L.tail headers) (L.tail body)]
                   "Data/plot_matrix"
                   "Plot Matrix"

    -- heatmap の作成はちょっと難しいのでpass

    -- 重回帰
    let (coeffs, r2) = olsRegress   [ VU.fromList index
                                    , VU.fromList hospital]
                                    ( VU.fromList death)

    putStrLn $ "回帰係数:" ++ show (L.init (VU.toList coeffs))
    putStrLn $ "切片:"     ++ show (L.last (VU.toList coeffs))

    -- 結果のプロット
    plotBar (L.zipWith  (\x y -> (x,[y]))
                        ["Economic Indicators", "Number of Hospitals"]
                        (L.init (VU.toList coeffs)))
            "Data/bar_result"
            "Plot Result"
            ["Coeff"]


text:
{-# LANGUAGE ForeignFunctionInterface #-}

-- $ ghc -O2 --make wrappers.hs functions.c

import Numeric.LinearAlgebra
import Data.Packed.Development
import Foreign(Ptr,unsafePerformIO)
import Foreign.C.Types(CInt)

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

main = do
    print $ myScale 3.0 (fromList [1..10])
    print $ myDiag $ (3><5) [1..]

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

foreign import ccall unsafe "c_scale_vector"
    cScaleVector :: Double                -- scale
                 -> CInt -> Ptr Double    -- argument
                 -> CInt -> Ptr Double    -- result
                 -> IO CInt               -- exit code

myScale s x = unsafePerformIO $ do
    y <- createVector (dim x)
    app2 (cScaleVector s) vec x vec y "cScaleVector"
    return y

-----------------------------------------------------
-- forcing row order

foreign import ccall unsafe "c_diag"
    cDiag :: CInt -> CInt -> Ptr Double  -- argument
          -> CInt -> Ptr Double          -- result1
          -> CInt -> CInt -> Ptr Double  -- result2
          -> IO CInt                     -- exit code

myDiag m = unsafePerformIO $ do
    y <- createVector (min r c)
    z <- createMatrix RowMajor r c
    app3 cDiag mat (cmat m) vec y mat z "cDiag"
    return (y,z)
  where r = rows m
        c = cols m

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
module FFT (fft) where
import Data.Complex (Complex, cis, conjugate)
import Data.List (head, zip, take)

fft :: [Complex Float] -> [Complex Float]
fft xs = fft' (length xs) xs

fft' :: Int -> [Complex Float] -> [Complex Float]
fft' 1 [x] = [x]
fft' n coeffs = eval [0..n-1] halfN nThPrimitiveRoots halfFFT
  where
        halfN                   = n `div` 2
        (evenCoeffs, oddCoeffs) = split coeffs
        evenFFT                 = fft' halfN evenCoeffs
        oddFFT                  = fft' halfN oddCoeffs
        nThPrimitiveRoots       = duplicate $ take halfN (rootsOfUnity n)
        halfFFT                 = duplicate $ zip evenFFT oddFFT

rootsOfUnity :: Int -> [Complex Float]
rootsOfUnity = rootsOfUnity' . fromIntegral

rootsOfUnity' :: Float -> [Complex Float]
rootsOfUnity' n = map conjugate $ [ cis (2 * pi * k / n) | k <- [0 .. n - 1] ]

split :: [Complex Float] -> ([Complex Float], [Complex Float])
split coeffs = (evenCoeffs, oddCoeffs)
  where indexedCoeffs = zip coeffs [0..]
        evenCoeffs    = [ coeff | (coeff, index) <- indexedCoeffs, even index ]
        oddCoeffs     = [ coeff | (coeff, index) <- indexedCoeffs, odd index ]

eval :: [Int] -> Int -> [Complex Float] -> [(Complex Float, Complex Float)] -> [Complex Float]
eval [] _ _ _ = []
eval (power:restPowers) halfPower (root:restRoots) (coeffs:restCoeffs)
  | power < halfPower = (evenCoeff + root * oddCoeff) : (eval restPowers halfPower restRoots restCoeffs)
  | otherwise         = (evenCoeff - root * oddCoeff) : (eval restPowers halfPower restRoots restCoeffs)
  where (evenCoeff, oddCoeff) = coeffs

duplicate :: [a] -> [a]
duplicate xs = xs ++ xs

max_stars_repo_name: gboduljak/dft

max_issues_repo_name: gboduljak/dft

max_forks_repo_name: gboduljak/dft


text:
{-# LANGUAGE OverloadedStrings #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_GHC -Wcompat #-}
{-# OPTIONS_GHC -Wincomplete-record-updates #-}
{-# OPTIONS_GHC -Wincomplete-uni-patterns #-}
{-# OPTIONS_GHC -Wredundant-constraints #-}

module CloudSunlight
  ( cloudSunlightExperiment,
  )
where

import Control.Monad
import qualified Data.ByteString.Lazy as BL
import Data.Csv
import qualified Data.List as L
import qualified Data.Vector.Unboxed as V
import qualified Statistics.Regression as R
import qualified Statistics.Sample as S
import System.Random.MWC
import qualified System.Random.MWC.Distributions as D

meanToa :: Double
meanToa = 340.0

toa :: GenIO -> IO Double
toa = D.normal meanToa 30.0

-- | magnitude of aerosol affect on sunlight, between 1.0 and 0.0 (1.0 = max effect)
aerosolAffectSunlight :: Double
aerosolAffectSunlight = 1.0

-- | factor that cloud reduces sunlight, if there is cloud
cloudAffectSunlight :: Double
cloudAffectSunlight = 0.6

sunlight ::
  -- | non-dim measure of aerosol
  Double ->
  -- | cloud
  Bool ->
  -- | U_sun (i.e. TOA)
  Double ->
  -- | W/m2
  Double
sunlight a c toa' = if c then clearSky * cloudAffectSunlight else clearSky
  where
    clearSky = toa' * (1.0 - aerosolAffectSunlight * a)

cloud :: Double -> Double -> Bool
cloud a u
  | (a + u) > 1.0 = True
  | otherwise = False

aerosol :: Double -> Double
aerosol u = u

genTuple :: GenIO -> IO (Double, String, Double)
genTuple gen = do
  uAero <- uniform gen
  uCld <- uniform gen
  uSun <- toa gen
  let a = aerosol uAero
  let c = cloud a uCld
  let s = sunlight a c uSun
  return (a, if c then "Cloudy" else "Clear", s)

getCloud :: (a, b, c) -> b
getCloud (_a, c, _s) = c

getSun :: (a, b, c) -> c
getSun (_a, _c, s) = s

getAerosol :: (a, b, c) -> a
getAerosol (a, _c, _s) = a

calcRegression :: V.Vector Double -> Double -> Double
calcRegression ols x = x * ols V.! 0 + ols V.! 1

calcStats ::
  [(Double, String, Double)] ->
  [(Double, String, Double)] ->
  (V.Vector Double, V.Vector Double, Double, Double)
calcStats cloud' clear =
  ( olsCloud,
    olsClear,
    avgRegress olsCloud - avgRegress olsClear,
    naiveCloud - naiveClear
  )
  where
    calcAvgs xs = (as, ols, S.mean ss)
      where
        ss = V.fromList (fmap getSun xs)
        as = V.fromList (fmap getAerosol xs)
        ols = fst . R.olsRegress [as] $ ss
    (aCloud, olsCloud, naiveCloud) = calcAvgs cloud'
    (aClear, olsClear, naiveClear) = calcAvgs clear
    totAs = aCloud <> aClear
    avgRegress ols = S.mean $ V.map (calcRegression ols) totAs

cloudSunlightExperiment :: Int -> IO (Double, Double, Double)
cloudSunlightExperiment n = do
  gen <- create
  xs <- replicateM n (genTuple gen)
  let (cloud', clear) = L.partition ((== "Cloudy") . getCloud) xs
  let (olsCloud, olsClear, estDiff, naiveDiff) = calcStats cloud' clear
  let regPoints = [0, 0.1 .. 1]
  f "dat/cloudy.dat" cloud'
  f "dat/clear.dat" clear
  write
    "aerosol\tsunlight\n"
    "dat/olsCloudy.dat"
    (zip regPoints (fmap (calcRegression olsCloud) regPoints))
  write
    "aerosol\tsunlight\n"
    "dat/olsClear.dat"
    (zip regPoints (fmap (calcRegression olsClear) regPoints))
  return (estDiff, naiveDiff, (0.5 * meanToa * (cloudAffectSunlight - 1.0)))
  where
    f = write "aerosol\tcloud\tsunlight\n"
    write header' fpath =
      BL.writeFile fpath . (header' <>)
        . encodeWith
          (defaultEncodeOptions {encDelimiter = 9})

max_stars_repo_name: massma/causality-earth-science

max_issues_repo_name: massma/causality-earth-science

max_forks_repo_name: massma/causality-earth-science


text:
-- | A basic interface between Vectors and the fftw library.
module Numeric.FFT.Vector.Base(
            -- * Transforms
            Transform(..),
            planOfType,
            PlanType(..),
            plan,
            run,
            -- * multi-demensional Transforms
            TransformND(..),
            planOfTypeND,
            planND,
            runND,
            -- * Plans
            Plan(),
            planInputSize,
            planOutputSize,
            execute,
            executeM,
            withPlanner,
            -- * Unsafe C stuff
            CFlags,
            CPlan,
            -- * Normalization helpers
            Scalable(..),
            modifyInput,
            modifyOutput,
            constMultOutput,
            multC,
            unsafeModify,
            ) where

import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as MS
import Data.Vector.Generic as V hiding (forM_)
import Data.Vector.Generic.Mutable as M hiding (unsafeModify)
import Data.List as L
import Control.Concurrent.MVar
import Control.Monad.Primitive (RealWorld,PrimMonad(..), PrimBase,
            unsafePrimToPrim, unsafePrimToIO)
import Control.Monad(forM_)
import Foreign (Storable(..), Ptr, FunPtr,
                ForeignPtr, withForeignPtr, newForeignPtr)
import Foreign.C (CInt(..), CUInt, CSize(..))
import Data.Bits ( (.|.) )
import Data.Complex(Complex(..))
import Foreign.Storable.Complex()
import System.IO.Unsafe (unsafePerformIO)

#include <fftw3.h>

---------------------
-- Creating FFTW plans

-- First, the Transform flags:
data PlanType = Estimate | Measure | Patient | Exhaustive
data Preservation = PreserveInput | DestroyInput

type CFlags = CUInt

-- | Marshal the Transform flags for use by fftw.
planInitFlags :: PlanType -> Preservation -> CFlags
planInitFlags pt pr = planTypeInt .|. preservationInt
  where
    planTypeInt = case pt of
                    Estimate -> #const FFTW_ESTIMATE
                    Measure -> #const FFTW_MEASURE
                    Patient -> #const FFTW_PATIENT
                    Exhaustive -> #const FFTW_EXHAUSTIVE
    preservationInt = case pr of
                    PreserveInput -> #const FFTW_PRESERVE_INPUT
                    DestroyInput -> #const FFTW_DESTROY_INPUT

newtype CPlan = CPlan {unCPlan :: ForeignPtr CPlan}

withPlan :: CPlan -> (Ptr CPlan -> IO a) -> IO a
withPlan = withForeignPtr . unCPlan

foreign import ccall unsafe fftw_execute :: Ptr CPlan -> IO ()
foreign import ccall "&" fftw_destroy_plan :: FunPtr (Ptr CPlan -> IO ())

newPlan :: Ptr CPlan -> IO CPlan
newPlan = fmap CPlan . newForeignPtr fftw_destroy_plan

----------------------------------------
-- vector-fftw plans

-- | A 'Plan' can be used to run an @fftw@ algorithm for a specific input/output size.
data Plan a b = Plan {
                    planInput :: {-# UNPACK #-} !(VS.MVector RealWorld a),
                    planOutput :: {-# UNPACK #-} !(VS.MVector RealWorld b),
                    planExecute :: IO ()
                }

-- | The (only) valid input size for this plan.
planInputSize :: Storable a => Plan a b -> Int
planInputSize = MS.length . planInput

-- | The (only) valid output size for this plan.
planOutputSize :: Storable b => Plan a b -> Int
planOutputSize = MS.length . planOutput

-- | Run a plan on the given 'Vector'.
--
-- If @'planInputSize' p /= length v@, then calling
-- @execute p v@ will throw an exception.
execute :: (Vector v a, Vector v b, Storable a, Storable b)
            => Plan a b -> v a -> v b
execute Plan{..} = \v -> -- fudge the arity to make sure it's always inlined
    if n /= V.length v
        then error $ "execute: size mismatch; expected " L.++ show n
                    L.++ ", got " L.++ show (V.length v)
        else unsafePerformIO $ do
                        forM_ [0..n-1] $ \k -> M.unsafeWrite planInput k
                                                $ V.unsafeIndex v k
                        planExecute
                        v' <- unsafeNew m
                        forM_ [0..m-1] $ \k -> M.unsafeRead planOutput k
                                                >>= M.unsafeWrite v' k
                        V.unsafeFreeze v'
  where
    n = MS.length planInput
    m = MS.length planOutput
{-# INLINE execute #-}

-- TODO: decide whether this is actually unsafe.
-- | Run a plan on the given mutable vectors.  The same vector may be used for both
-- input and output.
--
-- If @'planInputSize' p \/= length vIn@ or @'planOutputSize' p \/= length vOut@,
-- then calling @unsafeExecuteM p vIn vOut@ will throw an exception.
executeM :: forall m v a b .
        (PrimBase m, MVector v a, MVector v b, Storable a, Storable b)
            => Plan a b -- ^ The plan to run.
            -> v (PrimState m) a  -- ^ The input vector.
                    -> v (PrimState m) b -- ^ The output vector.
                    -> m ()
executeM Plan{..} = \vIn vOut ->
    if n /= M.length vIn || m /= M.length vOut
        then error $ "executeM: size mismatch; expected " L.++ show (n,m)
                    L.++ ", got " L.++ show (M.length vIn, M.length vOut)
        else unsafePrimToPrim $ act vIn vOut
  where
    n = MS.length planInput
    m = MS.length planOutput

    act :: v (PrimState m) a -> v (PrimState m) b -> IO ()
    act vIn vOut = do
            forM_ [0..n-1] $ \k -> unsafePrimToIO (M.unsafeRead vIn k :: m a)
                                    >>= M.unsafeWrite planInput k
            unsafePrimToPrim planExecute
            forM_ [0..n-1] $ \k -> M.unsafeRead planOutput k
                                    >>= unsafePrimToIO . (M.unsafeWrite vOut k
                                                            :: b -> m ())
{-# INLINE executeM #-}

------------------
-- Malloc/free of fftw array

foreign import ccall unsafe fftw_malloc :: CSize -> IO (Ptr a)
foreign import ccall "&" fftw_free :: FunPtr (Ptr a -> IO ())

newFFTVector :: forall a . Storable a => Int -> IO (MS.MVector RealWorld a)
newFFTVector n = do
    p <- fftw_malloc $ toEnum $ n * sizeOf (undefined :: a)
    fp <- newForeignPtr fftw_free p
    return $ MS.MVector n fp
{-# INLINE newFFTVector #-}

-----------------------
-- Transforms: methods of plan creation.

-- | A transform which may be applied to vectors of different sizes.
data Transform a b = Transform {
                        inputSize :: Int -> Int,
                        outputSize :: Int -> Int,
                        creationSizeFromInput :: Int -> Int,
                        makePlan :: CInt -> Ptr a -> Ptr b -> CFlags -> IO (Ptr CPlan),
                        normalization :: Int -> Plan a b -> Plan a b
                    }

-- | Create a 'Plan' of a specific size for this transform.
planOfType :: (Storable a, Storable b) => PlanType
                                -> Transform a b -> Int -> Plan a b
planOfType ptype Transform{..} n
  | m_in <= 0 || m_out <= 0 = error "Can't (yet) plan for empty arrays!"
  | otherwise  = unsafePerformIO $ do
    planInput <- newFFTVector m_in
    planOutput <- newFFTVector m_out
    MS.unsafeWith planInput $ \inP -> MS.unsafeWith planOutput $ \outP -> do
        pPlan <- makePlan (toEnum n) inP outP $ planInitFlags ptype DestroyInput
        cPlan <- newPlan pPlan
        -- Use unsafeWith here to ensure that the Storable MVectors' ForeignPtrs
        -- aren't released too soon:
        let planExecute = MS.unsafeWith planInput $ \_ ->
                            MS.unsafeWith planOutput $ \_ ->
                                withPlan cPlan fftw_execute
        return $ normalization n $ Plan {..}
  where
    m_in = inputSize n
    m_out = outputSize n
{-# INLINE planOfType #-}

-- | Create a 'Plan' of a specific size.  This function is equivalent to
-- @'planOfType' 'Estimate'@.
plan :: (Storable a, Storable b) => Transform a b -> Int -> Plan a b
plan = planOfType Estimate
{-# INLINE plan #-}

-- | Create and run a 'Plan' for the given transform.
run :: (Vector v a, Vector v b, Storable a, Storable b)
            => Transform a b -> v a -> v b
run p = \v -> execute
            (planOfType Estimate p $ creationSizeFromInput p $ V.length v)
            v
{-# INLINE run #-}

-----------------------
-- TransformND: methods of plan creation for multi-dimensional plans.

-- | A transform which may be applied to vectors of different sizes.
--
-- @since 0.2
data TransformND a b = TransformND {
                        inputSizeND :: Int -> Int,
                        outputSizeND :: Int -> Int,
                        creationSizeFromInputND :: Int -> Int,
                        makePlanND :: CInt -> Ptr CInt -> Ptr a -> Ptr b -> CFlags -> IO (Ptr CPlan),
                        normalizationND :: VS.Vector Int -> Plan a b -> Plan a b
                    }

-- | Create a 'Plan' of a specific size for this transform.
-- 'dims' must have rank greater or equal to 1
--
-- @since 0.2
planOfTypeND :: (Storable a, Storable b) => PlanType
                                -> TransformND a b -> VS.Vector Int -> Plan a b
planOfTypeND ptype TransformND{..} dims
  | m_in <= 0 || m_out <= 0 = error "Can't (yet) plan for empty arrays!"
  | otherwise  = unsafePerformIO $ do
    mdims <- unsafeThaw $ V.map toEnum dims
    planInput <- newFFTVector m_in
    planOutput <- newFFTVector m_out
    MS.unsafeWith mdims $ \dimsP -> MS.unsafeWith planInput $ \inP -> MS.unsafeWith planOutput $ \outP -> do
        pPlan <- makePlanND (toEnum $ V.length dims) dimsP inP outP $ planInitFlags ptype DestroyInput
        cPlan <- newPlan pPlan
        -- Use unsafeWith here to ensure that the Storable MVectors' ForeignPtrs
        -- aren't released too soon:
        let planExecute = MS.unsafeWith planInput $ \_ ->
                                MS.unsafeWith planOutput $ \_ ->
                                withPlan cPlan fftw_execute
        return $ normalizationND dims $ Plan {..}
  where
    m = V.product $ V.init dims
    m_in = m * inputSizeND (V.last dims)
    m_out = m * outputSizeND (V.last dims)
{-# INLINE planOfTypeND #-}

-- | Create a 'Plan' of a specific size.  This function is equivalent to
-- @'planOfType' 'Estimate'@.
--
-- @since 0.2
planND :: (Storable a, Storable b) => TransformND a b -> VS.Vector Int -> Plan a b
planND = planOfTypeND Estimate
{-# INLINE planND #-}

-- | Create and run a 'Plan' for the given transform.
--
-- @since 0.2
runND :: (Vector v a, Vector v b, Storable a, Storable b)
            => TransformND a b -> VS.Vector Int ->  v a -> v b
runND p = \dims v ->
  let creationSize = V.init dims `V.snoc` creationSizeFromInputND p (V.last dims) in
    execute
    (planOfTypeND Estimate p creationSize)
    v
{-# INLINE runND #-}

---------------------------
-- For scaling input/output:

class Scalable a where
    scaleByD :: Double -> a -> a

instance Scalable Double where
    scaleByD = (*)
    {-# INLINE scaleByD #-}

instance Scalable (Complex Double) where
    scaleByD s (x:+y) = s*x :+ s*y
    {-# INLINE scaleByD #-}

{-# INLINE modifyInput #-}
modifyInput :: (MS.MVector RealWorld a -> IO ()) -> Plan a b -> Plan a b
modifyInput f p@Plan{..} = p {planExecute = f planInput >> planExecute}

{-# INLINE modifyOutput #-}
modifyOutput :: (MS.MVector RealWorld b -> IO ()) -> Plan a b -> Plan a b
modifyOutput f p@Plan{..} = p {planExecute = planExecute >> f planOutput}

{-# INLINE constMultOutput #-}
constMultOutput :: (Storable b, Scalable b) => Double -> Plan a b -> Plan a b
constMultOutput !s = modifyOutput (multC s)

{-# INLINE multC #-}
multC :: (Storable a, Scalable a) => Double -> MS.MVector RealWorld a -> IO ()
multC !s v = forM_ [0..n-1] $ \k -> unsafeModify v k (scaleByD s)
  where !n = MS.length v

-- | Helper function; seems like it should be in the vector package...
{-# INLINE unsafeModify #-}
unsafeModify :: (Storable a)
                => MS.MVector RealWorld a -> Int -> (a -> a) -> IO ()
unsafeModify v k f = MS.unsafeRead v k >>= MS.unsafeWrite v k . f

plannerLock :: MVar ()
plannerLock = unsafePerformIO $ newMVar ()
{-# NOINLINE plannerLock #-}

-- | Calls to the FFTW planner are non-reentrant. Here we take a mutex to
-- ensure thread safety.
withPlanner :: IO a -> IO a
withPlanner action = do
    takeMVar plannerLock
    res <- action
    putMVar plannerLock ()
    return res

ext: hsc

max_stars_repo_name: TravisWhitaker/vector-fftw

max_issues_repo_name: TravisWhitaker/vector-fftw

max_forks_repo_name: TravisWhitaker/vector-fftw


text:
{-# LANGUAGE FlexibleContexts, NamedFieldPuns #-}

module School.Unit.CostFunction
( SetupCost
, CostFunction(..)
, defSetupCost
) where

import Conduit ((.|), ConduitM, ZipSource(..), getZipSource, mapC, mapMC)
import Numeric.LinearAlgebra (Container, Matrix, Vector, add, cols, rows)
import School.Types.Slinky (Slinky(..), slinkyConcat)
import School.Unit.CostParams (CostParams)
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitForward (ForwardStack)
import School.Unit.UnitGradient (UnitGradient(..))
import School.Utils.LinearAlgebra (zeroMatrix)

type SetupCost a m = ConduitM () (Matrix a) m ()
                  -> ConduitM () (ForwardStack a) m ()

defSetupCost :: (Monad m) => SetupCost a m
defSetupCost source =
  source .| mapMC (\matrix -> return ([BatchActivation matrix], SNil))

type CompupteCost a = UnitActivation a
                   -> Slinky CostParams
                   -> Either String a

type DerivCost a = UnitActivation a
                -> Slinky CostParams
                -> Either String (UnitGradient a)

type PrepareCost m a =
  ConduitM (ForwardStack a) (ForwardStack a) m ()

data CostFunction a m =
  CostFunction { computeCost :: CompupteCost a
               , derivCost :: DerivCost a
               , prepareCost :: PrepareCost m a
               , setupCost :: SetupCost a m
               }

mappendCompute :: (Num a)
               => CompupteCost a
               -> CompupteCost a
               -> CompupteCost a
mappendCompute c1 c2 activations params1@(SNode _ params2) = do
  res1 <- c1 activations params1
  res2 <- c2 activations params2
  return $ res1 + res2
mappendCompute _ _ _ _ = Left "mappend costfunctions needs two cost params"

mappendDeriv :: (Container Vector a, Num a)
             => DerivCost a
             -> DerivCost a
             -> DerivCost a
mappendDeriv d1 d2 activations params1@(SNode _ params2) = do
  (BatchGradient res1) <- d1 activations params1
  (BatchGradient res2) <- d2 activations params2
  return $ BatchGradient $ add res1 res2
mappendDeriv _ _ _ _ = Left "mappend costfunctions needs two cost params"

combinator :: ForwardStack a -> ForwardStack a -> ForwardStack a
combinator (_, cParams1) (activations, cParams2) =
  (activations, slinkyConcat cParams1 cParams2)

mappendSetup :: (Monad m)
             => SetupCost a m
             -> SetupCost a m
             -> SetupCost a m
mappendSetup s1 s2 source = getZipSource $ combinator
                                       <$> (ZipSource $ s1 source)
                                       <*> (ZipSource $ s2 source)

instance (Container Vector a, Num a, Monad m) => Monoid (CostFunction a m) where
  mappend CostFunction { computeCost = c1
                       , derivCost = d1
                       , prepareCost = p1
                       , setupCost = s1
                       }
          CostFunction { computeCost = c2
                       , derivCost = d2
                       , prepareCost = p2
                       , setupCost = s2
                       } =
   CostFunction { computeCost = mappendCompute c1 c2
                , derivCost = mappendDeriv d1 d2
                , prepareCost = p1 .| p2
                , setupCost = mappendSetup s1 s2
                }
  mempty = CostFunction { computeCost
                        , derivCost
                        , prepareCost = mapC id
                        , setupCost = defSetupCost
                        } where
    computeCost _ _ = Right 0
    derivCost (BatchActivation g) _ =
      Right . BatchGradient $ zeroMatrix r c
      where r = rows g
            c = cols g
    derivCost _ _ = Left "mempty deriv cost expects batch activation"

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
module Fractals
( order
, floatOrder
, mandelbrot
, bmp
, tangentPoint
, mandelbrotCardioid
, mag2
, toBitmap
) where

import Data.Complex (Complex(..), mkPolar)
import Graphics.Gloss
import qualified Data.ByteString as B
import Data.Word
import Control.Parallel.Strategies

import Palette

mag2:: Complex Double -> Double
mag2 (a:+b) = a^(2::Int) + b^(2::Int)

-- | Compute the order of a given complex function.
order :: (Complex Double -> Complex Double) -- ^ a function on the complex plane
      -> Double                             -- ^ the escape radius squared
      -> Int                                -- ^ the maximum order
      -> Complex Double                     -- ^ a starting point
      -> (Int, Complex Double)              -- ^ the order
order f radius2 nMax z0 = let f' = \ (i, z) -> (i+1, f z)
                              fs = dropWhile (\ (i, z) -> i<nMax && mag2 z < radius2) $ iterate f' (0, z0)
                           in head fs

-- | Compute the floating order of a given complex function.
floatOrder :: (Complex Double -> Complex Double) -- ^ a function on the complex plane
           -> Double                             -- ^ the escape radius squared
           -> Int                                -- ^ the maximum order
           -> Complex Double                     -- ^ a starting point
           -> Double                             -- ^ the order
floatOrder f radius2 nMax z0 = let (nu, z) = order f radius2 nMax z0
                                in fromIntegral nu - log (0.5 * (log $ mag2 z)) / log 2.0

-- | the Mandelbrot function: z^2 + c.
mandelbrot :: Complex Double -- ^ the c parameter
           -> Complex Double -- ^ the z variable
           -> Complex Double -- ^ the resulting value
mandelbrot (cr:+ci) (zr:+zi) = (cr + zr^(2::Int) - zi^(2::Int)) :+ (ci + 2*zr*zi)
--mandelbrot c z = z**2 + c

-- | Yield the point tangent to the main cardioid and to the bulb with period q
--   and combinatorial rotation number r=p/q.
tangentPoint :: Rational -> Complex Double
tangentPoint r = mandelbrotCardioid $ mkPolar 1.0 (2.0 * pi * fromRational r)

-- | The equation of the main cardioid.
mandelbrotCardioid :: Complex Double -> Complex Double
mandelbrotCardioid mu = 0.5 * mu * (1.0 - 0.5 * mu)

bitmapFormat :: BitmapFormat
bitmapFormat = BitmapFormat BottomToTop PxRGBA

escapeRadius2 :: Double
escapeRadius2 = 100000.0

maxIter :: Int
maxIter = 255

bmpValues :: Int -> Int -> Double -> Double -> Double -> [Double]
bmpValues w h xc yc picScale =
    let wf = fromIntegral w * picScale
        hf = fromIntegral h * picScale
        floatOrderFromIJ :: (Int, Int) -> Double
        floatOrderFromIJ (i,j) = let x0 = xc - 0.5*wf + (fromIntegral j)*wf/(fromIntegral (w-1))
                                     y0 = yc - 0.5*hf + (fromIntegral i)*hf/(fromIntegral (h-1))
                                     z0 = x0 :+ y0
                                  in floatOrder (mandelbrot z0) escapeRadius2 maxIter z0
        indices = (,) <$> [0..h-1] <*> [0..w-1]
        vals = map floatOrderFromIJ indices
     in vals `using` parListChunk h rseq

toBitmap :: Palette -> [Double] -> [Word8]
toBitmap _ [] = []
toBitmap palette l = let maxL = maximum l
                         toBitmap' _ [] = []
                         toBitmap' maxX (x:xs) = let (r,g,b) = palette (x/maxX)
                                                     in r:g:b:255:(toBitmap' maxX xs)
                      in toBitmap' maxL l

bmpByteString :: Int -> Int -> Double -> Double -> Double -> B.ByteString
bmpByteString w h xc yc picScale = B.pack vals
    where vals = toBitmap hsv $ bmpValues w h xc yc picScale

bmp :: Int -> Int -> Complex Double -> Double -> Picture
bmp w h (xc :+ yc) picScale = bitmapOfByteString w h bitmapFormat (bmpByteString w h xc yc picScale) False

max_stars_repo_name: arekfu/fractals-haskell

max_issues_repo_name: arekfu/fractals-haskell

max_forks_repo_name: arekfu/fractals-haskell


text:
module Arrow where

import           Graphics.Gloss
import           Data.Complex

type Arrow = Complex Float

vectorize :: Arrow -> Vector
vectorize a = (realPart a, imagPart a)

-- Rotate an arrow according to its frequency in cycles per second
updateArrow :: Arrow -> Integer -> Float -> Arrow
updateArrow arrow frequency t =
  cis (2 * pi * fromInteger frequency * t) * arrow

-- Create a path through a list of arrows arranged tip to tail
arrange :: [Arrow] -> Path
arrange = map vectorize . scanl (+) (0 :+ 0)

endpoint :: [Arrow] -> Point
endpoint = last . arrange

max_stars_repo_name: juliagracedefoor/fourier-visualizer

max_issues_repo_name: juliagracedefoor/fourier-visualizer

max_forks_repo_name: juliagracedefoor/fourier-visualizer


text:
module Tools
( randomMatrix
, makeBatchList
, encorderForDat2
, encorderForDat
, length'
, pickUpNum
, setElem
, oneHotEncoder
, oneHotEncoder'
) where

-------------------------------------------------------------
--module
-------------------------------------------------------------

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

import Durstenfeld

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

-------------------------------------------------------------
--Tools
-------------------------------------------------------------

--重みの初期化用
--要素が0。0から1.0の間の数で、指定した行数、列数のMatrixを生成する。
randomMatrix :: RandomGen g => Int -> Int -> g -> Matrix Double
randomMatrix row col g = matrix col $ take (row*col) $ randomRs (0.0, 0.01) g

--ex) randomMatrix 3 4 (mkStdGen 100)

makeBatchList :: (Eq t, Num t) => Int -> t -> Int -> Int -> [[Int]]
makeBatchList _ 0 _ _ = []
makeBatchList g num ds bs = batch' : makeBatchList (g - 1) (num - 1) ds bs
    where data_size = ds
          batch_size = bs
          batch' = randomChoice data_size batch_size (mkStdGen g )

--数字のリストのリストをファイルに書き込むために、良い感じにいじっている。
encorderForDat :: Show a => [[a]] -> [Char]
encorderForDat xs = intercalate "\n" $  fmap (intercalate " " . fmap show) xs

--リストを書き込む用
encorderForDat2 :: Show a => [a] -> [Char]
encorderForDat2 xs = intercalate "\n" $  fmap show xs

length' :: [a] -> Double
length' [] = 0.0
length' (x:xs) = 1.0 + length' xs

pickUpNum :: [Int] -> [a] -> [a]
pickUpNum [] _ = []
pickUpNum (x:xs) list = number : pickUpNum xs list
    where number = list !! x

setElem1 :: Matrix R -> (Int, Int) -> Int -> Int -> R -> Matrix R
setElem1 m (row, col) maxrow maxcol newElem = newM
    where v  = flatten m
          num = row * maxcol + (col+1) --要素が何番目か
          v1 = subVector 0 (num-1) v
          v2 = subVector (num) (maxcol*maxrow-num) v
          v3 = vector [newElem]
          newV = vjoin [v1,v3,v2]
          newM = reshape maxcol newV

--〇
setElem :: Matrix R -> (Int, Int) -> R -> Matrix R
setElem m (row, col) newElem = setElem1 m (row, col) maxrow maxcol newElem
    where maxrow = rows m
          maxcol = cols m

--ラベル(0~9)をOne-Hot表現のMatrix(1行10列)に変換
--ex) 9 --> 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
oneHotEncoder :: Int -> Matrix Double
oneHotEncoder num = setElem zeroM (0 ,num) 1.0
    where zeroM = row [0.0,0,0,0,0,0,0,0,0,0]

-- =============================================

--setElem1' :: Matrix R -> (Int, Int) -> Int -> Int -> R -> Matrix R
setElem' v col newElem = newV
    where v1 = subVector 0 (col-1) v
          v2 = subVector col (10-col) v
          v3 = vector [newElem]
          newV = vjoin [v1,v3,v2]

oneHotEncoder' :: Int -> Vector R
oneHotEncoder' num = setElem' zeroV (num+1) 1.0
    where zeroV = fromList [0.0,0,0,0,0,0,0,0,0,0]

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

max_stars_repo_name: llbxg/Fukami

max_issues_repo_name: llbxg/Fukami

max_forks_repo_name: llbxg/Fukami


text:
{-# LANGUAGE Strict, StrictData #-}
{-# LANGUAGE DeriveGeneric, DeriveAnyClass #-}
{-# LANGUAGE TypeFamilies, TypeFamilyDependencies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE DefaultSignatures #-}

------------------------------------------------
-- |
-- Module    :  Numeric.Algebra.Presentation
-- Copyright :  (c) Jun Yoshida 2019
-- License   :  BSD3
--
-- Matrix presentations of morphisms
--
------------------------------------------------

module Numeric.Algebra.Presentation where

import GHC.Generics (Generic)
import Control.DeepSeq (NFData, force)

import Control.Monad
import Control.Monad.ST

import Data.Foldable

import qualified Numeric.LinearAlgebra as LA
import qualified Numeric.LinearAlgebra.Devel as LA

import Numeric.Algebra.FreeModule
import Numeric.Matrix.Integral
import Numeric.F2 (F2)
import qualified Numeric.Matrix.F2 as F2
import qualified Numeric.Matrix.F2.Mutable as F2

---------------
-- * Classes
---------------
-- | Class which provides matrix presentations for homomorphisms with coefficients in a.
class (Eq a, Num a, NFData a, Semigroup (Mat a), NFData (Vec a), NFData (Mat a)) => Coefficient a where
  type Vec a = v | v -> a
  -- ^ The type of vectors with coefficients in a.
  type Mat a = m | m -> a
  -- ^ The type of matrices with coefficients in a.
  ident :: Int -> Mat a
  -- ^ Presentation of the identity of the given rank.
  present :: (Integral a', NFData a', Ord c, NFData c) => (b->FreeMod a' c) -> [b] -> [c] -> Mat a
  -- ^ Present a linear map in a matrix form with respect to a given bases
  vecToList :: Vec a -> [a]
  -- ^ Vec a is convertible to a list.
  rankVec :: Vec a -> Int
  -- ^ The rank of the module where the given vector lies.
  rankDom :: Mat a -> Int
  -- ^ The rank of the domain of the given matrix.
  rankCod :: Mat a -> Int
  -- ^ The rank of the codomain of the given matrix.
  fromVecs :: (Foldable t) => t (Vec a) -> Mat a
  -- ^ Give a matrix whose image is spanned by given vectors.
  -- Typically, it is the matrix with given vectors as columns.
  toVecs :: Mat a -> [Vec a]
  -- ^ Give a set of vectors which span the image of the matrix.
  -- Typically, it is the set of column vectors.
  -- ^ The scalar matrix of the given size.
  (@|>) :: Mat a -> Vec a -> Vec a
  -- ^ Multiplication of matrices and vectors.
  scale :: a -> Vec a -> Vec a
  -- ^ Scalar multiplication on vectors.
  scalar :: Int -> a -> Mat a
  -- ^ Scalar * identity matrix.
  -- | Default implementations
  ident n = scalar n 1
  rankVec = length . vecToList
  (@|>) x = head . toVecs . (x<>) . fromVecs . (:[])
  scale x v = scalar (rankVec v) x @|> v
  scalar n x = fromVecs (map (scale x) (toVecs (ident n)))

-- | Matrix presentations with normal forms.
class (Coefficient a, NFData (NormalForm a)) => Normalized a where
  data NormalForm a
  -- ^ The type of normal forms.
  fromNF :: NormalForm a -> Mat a
  -- ^ Convert into a matrix form.
  -- >> let xnf = toNF x
  -- Then
  -- >> x <> basisDom xnf == basisCod xnf <> fromNF xnf
  basisDom :: NormalForm a -> Mat a
  -- ^ The basis matrix on the domain.
  -- It should be invertible.
  basisCod :: NormalForm a -> Mat a
  -- ^ The basis matrix on the codomain.
  -- It should be invertible.
  toNF :: Mat a -> NormalForm a
  -- ^ Compute the normal form.
  invariants :: NormalForm a -> [a]
  -- ^ Each normal form should contain the data of "invariant factors;" e.g. those for Smith normal forms of matrices with coefficients in a PID.

------------------------------------------------
-- * Instances in terms of LA.Z using hmatrix
------------------------------------------------
instance Coefficient (LA.Z) where
  type Vec LA.Z = LA.Vector LA.Z
  type Mat LA.Z = LA.Matrix LA.Z
  ident = LA.ident
  present f dom cod
    = let domRk = length dom
          codRk = length cod
      in runST $ do
    stMat <- LA.newUndefinedMatrix LA.RowMajor codRk domRk
    forM_ [0..domRk-1] $ \j -> do
      let !vals = force $ f (dom !!j)
      forM_ [0..codRk-1] $ \i -> do
        let coeff = getCoeff (cod !! i) vals
        LA.unsafeWriteMatrix stMat i j (fromIntegral coeff)
    LA.unsafeFreezeMatrix stMat
  vecToList = LA.toList
  rankVec = LA.size
  rankDom = LA.cols
  rankCod = LA.rows
  fromVecs = LA.fromColumns . toList
  toVecs = LA.toColumns
  (@|>) = (LA.#>)
  scale = LA.scale

instance Normalized (LA.Z) where
  -- | Presentation in Smith normal forms.
  data NormalForm LA.Z = SmithNF {
    _invFactor :: LA.Vector LA.Z, -- Invariant factor
    _basisS :: LA.Matrix LA.Z,    -- basis for source
    _basisT :: LA.Matrix LA.Z     -- basis for target
    } deriving (Show, Generic, NFData)
  fromNF (SmithNF invs bss bts)
    = runST $ do
    stMat <- LA.newMatrix 0 (LA.cols bts) (LA.cols bss)
    forM_ [0..(LA.size invs)-1] $ \i ->
      LA.writeMatrix stMat i i  (invs `LA.atIndex` i)
    LA.unsafeFreezeMatrix stMat
  basisDom = _basisS
  basisCod = _basisT
  toNF x
    = let (u,s,v) = smithRep x
      in SmithNF (LA.takeDiag s) v u
  invariants = LA.toList . _invFactor

---------------------------------
-- * Instances in terms of F2
---------------------------------
instance Coefficient F2 where
  type Vec F2 = F2.VectorF2
  type Mat F2 = F2.MatrixF2
  ident = F2.ident
  present f dom cod
    = let domRk = length dom
          codRk = length cod
      in runST $ do
    stMat <- F2.unsafeNewMatrixF2 codRk domRk
    forM_ [0..domRk-1] $ \j -> do
      let !vals = force $ f (dom !!j)
      forM_ [0..codRk-1] $ \i -> do
        let coeff = getCoeff (cod !! i) vals
        F2.unsafeWriteMatrixF2 stMat i j (fromIntegral coeff :: F2)
    F2.unsafeFreezeMatrixF2 stMat
  vecToList = F2.toList
  rankVec = F2.sizeF2
  rankDom = F2.cols
  rankCod = F2.rows
  fromVecs = F2.fromColumns . toList
  toVecs = F2.toColumns
  scalar = F2.scalar

instance Normalized F2 where
  data NormalForm F2 =
    DiagNF {-# UNPACK #-} !Int {-# UNPACK #-} !F2.MatrixF2 {-# UNPACK #-} !F2.MatrixF2
    deriving (Show, Generic,NFData)
  invariants (DiagNF rk _ _) = replicate rk 1
  basisDom (DiagNF _ bss _) = bss
  basisCod (DiagNF _ _ bst) = bst
  fromNF (DiagNF rk bss bst) = F2.diagRectL (rankCod bst) (rankCod bss) (repeat 1 :: [Int])
  toNF x
    = let (u,rk,v) = F2.diagRep x
      in DiagNF rk v u

max_stars_repo_name: Junology/linkedgram

max_issues_repo_name: Junology/linkedgram

max_forks_repo_name: Junology/linkedgram