text:
{-# LANGUAGE CPP                                      #-}
{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE DeriveFoldable                           #-}
{-# LANGUAGE DeriveFunctor                            #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE FlexibleInstances                        #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE LambdaCase                               #-}
{-# LANGUAGE OverloadedStrings                        #-}
{-# LANGUAGE PatternSynonyms                          #-}
{-# LANGUAGE RecordWildCards                          #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE StandaloneDeriving                       #-}
{-# LANGUAGE TupleSections                            #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE TypeSynonymInstances                     #-}
{-# LANGUAGE ViewPatterns                             #-}
{-# OPTIONS_GHC -fno-warn-orphans                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

-- | Hamilton example suite
--
-- See: https://github.com/mstksg/hamilton#example-app-runner
--
-- Or just run with:
--
-- > $ hamilton-examples --help
-- > $ hamilton-examples [EXAMPLE] --help
--

import           Control.Concurrent
import           Control.Monad
import           Data.Bifunctor
import           Data.Finite
import           Data.Foldable
import           Data.IORef
import           Data.List
import           Data.Semigroup                      ((<>))
import           GHC.TypeLits
import           Graphics.Vty hiding                 (Config, (<|>))
import           Numeric.Hamilton
import           Numeric.LinearAlgebra.Static hiding (dim, (<>))
import           Numeric.LinearAlgebra.Static.Vector
import           Options.Applicative
import           System.Exit
import           Text.Printf
import           Text.Read
import qualified Data.List.NonEmpty                  as NE
import qualified Data.Map.Strict                     as M
import qualified Data.Vector                         as VV
import qualified Data.Vector.Sized                   as V
import qualified Data.Vector.Storable.Sized          as VS
import qualified Text.PrettyPrint.ANSI.Leijen        as PP

data SysExample where
    SE :: (KnownNat m, KnownNat n)
       => { seName   :: String
          , seCoords :: V.Vector n String
          , seSystem :: System m n
          , seDraw   :: R m -> [V2 Double]
          , seInit   :: Phase n
          }
       -> SysExample

pendulum :: Double -> Double -> SysExample
pendulum θ0 ω0 = SE "Single pendulum" (V1 "θ") s f (toPhase s c0)
  where
    s :: System 2 1
    s = mkSystem' (vec2 1 1                             )     -- masses
                  (\(V1 θ)   -> V2 (sin θ) (0.5 - cos θ))     -- coordinates
                  (\(V2 _ y) -> y                       )     -- potential
    f :: R 2 -> [V2 Double]
    f xs = [grVec xs]
    c0 :: Config 1
    c0 = Cfg (konst θ0 :: R 1) (konst ω0 :: R 1)

doublePendulum :: Double -> Double -> SysExample
doublePendulum m1 m2 = SE "Double pendulum" (V2 "θ1" "θ2") s f (toPhase s c0)
  where
    s :: System 4 2
    s = mkSystem' (vec4 m1 m1 m2 m2)     -- masses
                  (\(V2 θ1 θ2)     -> V4 (sin θ1)            (1 - cos θ1)
                                         (sin θ1 + sin θ2/2) (1 - cos θ1 - cos θ2/2)
                  )                      -- coordinates
                  (\(V4 _ y1 _ y2) -> 5 * (realToFrac m1 * y1 + realToFrac m2 * y2))
                                         -- potential
    f :: R 4 -> [V2 Double]
    f (split->(xs,ys))= grVec <$> [xs, ys]
    c0 :: Config 2
    c0 = Cfg (vec2 (pi/2) 0) (vec2 0 0)

room :: Double -> SysExample
room θ = SE "Room" (V2 "x" "y") s f (toPhase s c0)
  where
    s :: System 2 2
    s = mkSystem (vec2 1 1)         -- masses
                 id                 -- coordinates
                 (\(V2 x y) -> sum [ 2 * y                      -- gravity
                                   , 1 - logistic (-1) 10 0.1 y  -- bottom wall
                                   ,     logistic 1 10 0.1 y     -- top wall
                                   , 1 - logistic (-2) 10 0.1 x  -- left wall
                                   ,     logistic 2 10 0.1 x     -- right wall
                                   ]
                 )                  -- potential
    f :: R 2 -> [V2 Double]
    f xs = [grVec xs]
    c0 :: Config 2
    c0 = Cfg (vec2 (-1) 0.25) (vec2 (cos θ) (sin θ))

twoBody :: Double -> Double -> Double -> SysExample
twoBody m1 m2 ω0 = SE "Two-Body" (V2 "r" "θ") s f (toPhase s c0)
  where
    mT :: Double
    mT = m1 + m2
    s :: System 4 2
    s = mkSystem (vec4 m1 m1 m2 m2) -- masses
                 -- positions are calculated assuming (0,0) is the center
                 -- of mass
                 (\(V2 r θ) -> let r1 = r * realToFrac (-m2 / mT)
                                   r2 = r * realToFrac (m1 / mT)
                               in  V4 (r1 * cos θ) (r1 * sin θ)
                                      (r2 * cos θ) (r2 * sin θ)
                 )                 -- coordinates
                 (\(V2 r _) -> - realToFrac (m1 * m2) / r)  -- potential
    f :: R 4 -> [V2 Double]
    f (split->(xs,ys))= grVec <$> [xs, ys]
    c0 :: Config 2
    c0 = Cfg (vec2 2 0) (vec2 0 ω0)

spring
    :: Double -> Double -> Double -> Double -> SysExample
spring mB mW k x0 = SE "Spring hanging from block" (V3 "r" "x" "θ") s f (toPhase s c0)
  where
    s :: System 3 3
    s = mkSystem (vec3 mB mW mW)                                                  -- masses
                 (\(V3 r x θ)  -> V3 r (r + (1 + x) * sin θ) ((1 + x) * (-cos θ))) -- coordinates
                 (\(V3 r x θ) -> realToFrac k * x**2 / 2        -- spring
                              + (1 - logistic (-1.5) 25 0.1 r)  -- left rail wall
                              + (    logistic   1.5  25 0.1 r)  -- right rail wall
                              + realToFrac mB * ((1 + x) * (-cos θ))  -- gravity
                 )
    f :: R 3 -> [V2 Double]
    f (headTail->(b,w)) = [V2 b 1, V2 0 1 + grVec w]
    c0 :: Config 3
    c0 = Cfg (vec3 0 x0 0) (vec3 1 0 (-0.5))

bezier
    :: forall n. KnownNat (1 + n)
    => V.Vector (1 + n) (V2 Double)
    -> SysExample
bezier ps = SE "Bezier" (V1 "t") s f (toPhase s c0)
  where
    s :: System 2 1
    s = mkSystem (vec2 1 1)                                             -- masses
                 (\(V1 t) -> bezierCurve (fmap realToFrac <$> ps) t)    -- coordinates
                 (\(V1 t) -> (1 - logistic 0 5 0.05 t)           -- left wall
                           +      logistic 1 5 0.05 t            -- right wall
                 )
    f :: R 2 -> [V2 Double]
    f xs = [grVec xs]
    c0 :: Config 1
    c0 = Cfg (0.5 :: R 1) (0.25 :: R 1)

data ExampleOpts = EO { eoChoice :: SysExampleChoice }

data SysExampleChoice =
        SECDoublePend Double Double
      | SECPend Double Double
      | SECRoom Double
      | SECTwoBody Double Double Double
      | SECSpring Double Double Double Double
      | SECBezier (NE.NonEmpty (V2 Double))

parseEO :: Parser ExampleOpts
parseEO = EO <$> (parseSEC <|> pure (SECDoublePend 1 1))

parseSEC :: Parser SysExampleChoice
parseSEC = subparser . mconcat $
    [ command "doublepend" $
        info (helper <*> parseDoublePend)
             (progDesc "Double pendulum (default)")
    , command "pend"       $
        info (helper <*> parsePend      )
             (progDesc "Single pendulum")
    , command "room"       $
        info (helper <*> parseRoom      )
        (progDesc "Ball in room, bouncing off of walls")
    , command "twobody"    $
        info (helper <*> parseTwoBody    )
        (progDesc "Two-body graviational simulation.  Note that bodies will only orbit if H < 0.")
    , command "spring"    $
        info (helper <*> parseSpring    )
        (progDesc "A spring hanging from a block on a rail, holding up a mass.  Block is constrained to bounce between -1.5 and 1.5.")
    , command "bezier"     $
        info (helper <*> parseBezier    )
        (progDesc "Particle moving along a parameterized bezier curve")
    , metavar "EXAMPLE"
    ]
  where
    parsePend
      = SECPend       <$> option auto ( long "angle"
                                     <> short 'a'
                                     <> metavar "ANGLE"
                                     <> help "Intitial rightward angle (in degrees) of bob"
                                     <> value 0
                                     <> showDefault
                                      )
                      <*> option auto ( long "vel"
                                     <> short 'v'
                                     <> metavar "VELOCITY"
                                     <> help "Initial rightward angular velocity of bob"
                                     <> value 1
                                     <> showDefault
                                      )
    parseDoublePend
      = SECDoublePend <$> option auto ( long "m1"
                                     <> metavar "MASS"
                                     <> help "Mass of first bob"
                                     <> value 1
                                     <> showDefault
                                      )
                      <*> option auto ( long "m2"
                                     <> metavar "MASS"
                                     <> help "Mass of second bob"
                                     <> value 1
                                     <> showDefault
                                      )
    parseRoom
      = SECRoom    <$> option auto ( long "angle"
                                  <> short 'a'
                                  <> metavar "ANGLE"
                                  <> help "Initial upward launch angle (in degrees) of object"
                                  <> value 45
                                  <> showDefault
                                   )
    parseTwoBody
      = SECTwoBody <$> option auto ( long "m1"
                                  <> metavar "MASS"
                                  <> help "Mass of first body"
                                  <> value 5
                                  <> showDefault
                                   )
                   <*> option auto ( long "m2"
                                  <> metavar "MASS"
                                  <> help "Mass of second body"
                                  <> value 0.5
                                  <> showDefault
                                   )
                   <*> option auto ( long "vel"
                                  <> short 'v'
                                  <> metavar "VELOCITY"
                                  <> help "Initial angular velocity of system"
                                  <> value 0.5
                                  <> showDefault
                                   )
    parseSpring
      = SECSpring <$> option auto ( long "block"
                                 <> short 'b'
                                 <> metavar "MASS"
                                 <> help "Mass of block on rail"
                                 <> value 2
                                 <> showDefault
                                  )
                  <*> option auto ( long "weight"
                                 <> short 'w'
                                 <> metavar "MASS"
                                 <> help "Mass of weight hanging from spring"
                                 <> value 1
                                 <> showDefault
                                  )
                  <*> option auto ( short 'k'
                                 <> metavar "NUM"
                                 <> help "Spring constant / stiffness of spring"
                                 <> value 10
                                 <> showDefault
                                  )
                  <*> option auto ( short 'x'
                                 <> metavar "DIST"
                                 <> help "Initial displacement of spring"
                                 <> value 0.1
                                 <> showDefault
                                  )
    parseBezier
      = SECBezier <$> option f ( long "points"
                              <> short 'p'
                              <> metavar "POINTS"
                              <> help "List of control points (at least one), as tuples"
                              <> value (V2 (-1) (-1) NE.:| [V2 (-2) 1, V2 0 1, V2 1 (-1), V2 2 1])
                              <> showDefaultWith (show . map (\(V2 x y) -> (x, y)) . toList)
                               )
      where f = eitherReader $ \s -> do
              ps  <- maybe (Left "Bad parse") Right
                  $ readMaybe s
              maybe (Left "At least one control point required") Right
                  $ NE.nonEmpty (uncurry V2 <$> ps)

data SimOpts = SO { soZoom :: Double
                  , soRate :: Double
                  , soHist :: Int
                  }
  deriving (Show)

data SimEvt = SEQuit
            | SEZoom Double
            | SERate Double
            | SEHist Int

main :: IO ()
main = do
    EO{..} <- execParser $ info (helper <*> parseEO)
        ( fullDesc
       <> header "hamilton-examples - hamilton library example suite"
       <> progDescDoc (Just descr)
        )

    vty <- mkVty =<< standardIOConfig

    opts <- newIORef $ SO 0.5 1 25

    t <- forkIO . loop vty opts $ case eoChoice of
      SECDoublePend m1 m2        -> doublePendulum m1 m2
      SECPend       d0 ω0        -> pendulum (d0 / 180 * pi) ω0
      SECRoom       d0           -> room (d0 / 180 * pi)
      SECTwoBody    m1 m2 ω0     -> twoBody m1 m2 ω0
      SECSpring     mB mW k x0   -> spring mB mW k x0
      SECBezier     (p NE.:| ps) -> V.withSized (VV.fromList ps)
                                      (bezier . V.cons p)

    forever $ do
      e <- nextEvent vty
      forM_ (processEvt e) $ \case
        SEQuit -> do
          killThread t
          shutdown vty
          exitSuccess
        SEZoom s ->
          modifyIORef opts $ \o -> o { soZoom = soZoom o * s }
        SERate r ->
          modifyIORef opts $ \o -> o { soRate = soRate o * r }
        SEHist h ->
          modifyIORef opts $ \o -> o { soHist = soHist o + h }
  where
    fps :: Double
    fps = 12
    screenRatio :: Double
    screenRatio = 2.1
    ptAttrs :: [(Char, Color)]
    ptAttrs  = ptChars `zip` ptColors
      where
        ptColors = cycle [white,yellow,blue,red,green]
        ptChars  = cycle "o*+~"
    loop :: Vty -> IORef SimOpts -> SysExample -> IO ()
    loop vty oRef SE{..} = go M.empty seInit
      where
        qVec = intercalate "," . V.toList $ seCoords
        go hists p = do
          SO{..} <- readIORef oRef
          let p'   = stepHam (soRate / fps) seSystem p  -- progress the simulation
              xb   = (- recip soZoom, recip soZoom)
              infobox = vertCat . map (string defAttr) $
                  [ printf "[ %s ]" seName
                  , printf " <%s>   : <%s>" qVec . intercalate ", "
                     . map (printf "%.4f") . VS.toList . rVec . phsPositions $ p
                  , printf "d<%s>/dt: <%s>" qVec . intercalate ", "
                     . map (printf "%.4f") . VS.toList . rVec . velocities seSystem $ p
                  , printf "KE: %.4f" . keP seSystem           $ p
                  , printf "PE: %.4f" . pe seSystem . phsPositions $ p
                  , printf "H : %.4f" . hamiltonian seSystem   $ p
                  , " "
                  , printf "rate: x%.2f <>" $ soRate
                  , printf "hist: % 5d []" $ soHist
                  , printf "zoom: x%.2f -+" $ soZoom
                  ]
              pts  = (`zip` ptAttrs) . seDraw . underlyingPos seSystem . phsPositions
                   $ p
              hists' = foldl' (\h (r, a) -> M.insertWith (addHist soHist) a [r] h) hists pts
          dr <- displayBounds $ outputIface vty
          update vty . picForLayers . (infobox:) . plot dr (PX xb (RR 0.5 screenRatio)) $
               ((second . second) (defAttr `withForeColor`) <$> pts)
            ++ (map (\((_,c),r) -> (r, ('.', defAttr `withForeColor` c)))
                  . concatMap sequence
                  . M.toList
                  $ hists'
               )
          threadDelay (round (1000000 / fps))
          go hists' p'
    addHist hl new old = take hl (new ++ old)
    descr :: PP.Doc
    descr = PP.vcat
      [ "Run examples from the hamilton library example suite."
      , "Use with [EXAMPLE] --help for more per-example options."
      , ""
      , "To adjust rate/history/zoom, use keys <>/[]/-+, respectively."
      , ""
      , "See: https://github.com/mstksg/hamilton#example-app-runner"
      ]

processEvt
    :: Event -> Maybe SimEvt
processEvt = \case
    EvKey KEsc        []      -> Just SEQuit
    EvKey (KChar 'c') [MCtrl] -> Just SEQuit
    EvKey (KChar 'q') []      -> Just SEQuit
    EvKey (KChar '+') []      -> Just $ SEZoom (sqrt 2)
    EvKey (KChar '-') []      -> Just $ SEZoom (sqrt 0.5)
    EvKey (KChar '>') []      -> Just $ SERate (sqrt 2)
    EvKey (KChar '<') []      -> Just $ SERate (sqrt (1/2))
    EvKey (KChar ']') []      -> Just $ SEHist 5
    EvKey (KChar '[') []      -> Just $ SEHist (-5)
    _                         -> Nothing

data RangeRatio = RR { -- | Where on the screen (0 to 1) to place the other axis
                       rrZero  :: Double
                       -- | Ratio of height of a terminal character to width
                     , rrRatio :: Double
                     }
                deriving (Show)

data PlotRange = PXY (Double, Double) (Double, Double)
               | PX  (Double, Double) RangeRatio
               | PY  RangeRatio       (Double, Double)

plot
    :: (Int, Int)               -- ^ display bounds
    -> PlotRange
    -> [(V2 Double, (Char, Attr))]   -- ^ points to plot
    -> [Image]
plot (wd,ht) pr = map (crop wd ht)
                . (++ bgs)
                . map (\(p, (c, a)) -> place EQ EQ p $ char a c)
  where
    wd' = fromIntegral wd
    ht' = fromIntegral ht
    ((xmin, xmax), (ymin, ymax)) = mkRange (wd', ht') pr
    origin = place EQ EQ (V2 0 0) $ char defAttr '+'
    xaxis  = place EQ EQ (V2 0 0) $ charFill defAttr '-' wd 1
    yaxis  = place EQ EQ (V2 0 0) $ charFill defAttr '|' 1 ht
    xrange = xmax - xmin
    yrange = ymax - ymin
    bg     = backgroundFill wd ht
    scale (V2 pX pY) = V2 x y
      where
        x = round $ (pX - xmin) * (wd' / xrange)
        y = round $ (pY - ymin) * (ht' / yrange)
    place aX aY p i = case scale p of
      V2 pX pY -> translate (fAlign aX (imageWidth  i))
                            (fAlign aY (imageHeight i))
                . translate pX pY
                $ i
    labels = [ place LT EQ (V2 xmin 0) . string defAttr $ printf "%.2f" xmin
             , place GT EQ (V2 xmax 0) . string defAttr $ printf "%.2f" xmax
             , place EQ LT (V2 0 ymin) . string defAttr $ printf "%.2f" ymin
             , place EQ GT (V2 0 ymax) . string defAttr $ printf "%.2f" ymax
             ]
    bgs    = labels ++ [origin, xaxis, yaxis, bg]
    fAlign = \case
      LT -> const 0
      EQ -> negate . (`div` 2)
      GT -> negate

mkRange
    :: (Double, Double)
    -> PlotRange
    -> ((Double, Double), (Double, Double))
mkRange (wd, ht) = \case
    PXY xb     yb     -> (xb, yb)
    PX  xb     RR{..} ->
      let yr = (uncurry (-) xb) * ht / wd * rrRatio
          y0 = (rrZero - 1) * yr
      in  (xb, (y0, y0 + yr))
    PY  RR{..} yb ->
      let xr = (uncurry (-) yb) * wd / ht / rrRatio
          x0 = (rrZero - 1) * xr
      in  ((x0, x0 + xr), yb)

pattern V1 :: a -> V.Vector 1 a
pattern V1 x <- (V.head->x)
  where
    V1 x = V.singleton x
#if __GLASGOW_HASKELL__ >= 802
{-# COMPLETE V1 #-}
#endif

type V2 = V.Vector 2
pattern V2 :: a -> a -> V2 a
pattern V2 x y <- (V.toList->[x,y])
  where
    V2 x y = V.fromTuple (x, y)
#if __GLASGOW_HASKELL__ >= 802
{-# COMPLETE V2 #-}
#endif

pattern V3 :: a -> a -> a -> V.Vector 3 a
pattern V3 x y z <- (V.toList->[x,y,z])
  where
    V3 x y z = V.fromTuple (x, y, z)
#if __GLASGOW_HASKELL__ >= 802
{-# COMPLETE V3 #-}
#endif

pattern V4 :: a -> a -> a -> a -> V.Vector 4 a
pattern V4 x y z a <- (V.toList->[x,y,z,a])
  where
    V4 x y z a = V.fromTuple (x, y, z, a)
#if __GLASGOW_HASKELL__ >= 802
{-# COMPLETE V4 #-}
#endif

logistic
    :: Floating a => a -> a -> a -> a -> a
logistic pos ht width = \x -> ht / (1 + exp (- beta * (x - pos)))
  where
    beta = log (0.9 / (1 - 0.9)) / width

bezierCurve
    :: forall n f a. (KnownNat (1 + n), Applicative f, Num a)
    => V.Vector (1 + n) (f a)
    -> a
    -> f a
bezierCurve ps t =
      foldl' (liftA2 (+)) (pure 0)
    . V.imap (\i -> let i' = fromIntegral i
                    in  fmap (* (fromIntegral (n' `choose` i') * (1 - t)^(n' - i') * t^i))
             )
    $ ps
  where
    n' :: Int
    n' = fromIntegral (maxBound :: Finite (1 + n))
    choose :: Int -> Int -> Int
    n `choose` k = factorial n `div` (factorial (n - k) * factorial k)
    factorial :: Int -> Int
    factorial m = product [1..m]

deriving instance Ord Color

max_stars_repo_name: mstksg/hamilton

max_issues_repo_name: mstksg/hamilton

max_forks_repo_name: mstksg/hamilton


text:
module Scheme.Primitives.Numbers
  ( numericPrimitives,
    unwrapNumber,
  )
where

import Control.Exception (throw)
import Control.Monad.Except (MonadError (throwError))
import Data.Complex (imagPart, realPart)
import Data.Ratio (denominator, numerator)
import Data.Text (Text)
import Scheme.Operators (unOp)
import Scheme.Types (Number (..), SchemeError (..), SchemeResult, SchemeVal (..))

numericPrimitives :: [(Text, [SchemeVal] -> SchemeResult SchemeVal)]
numericPrimitives =
  [ ("number?", unOp isNumber),
    ("complex?", unOp isComplex),
    ("real?", unOp isReal),
    ("rational?", unOp isRational),
    ("integer?", unOp isInteger),
    ("exact?", unOp isExact),
    ("inexact?", unOp isInexact),
    ("exact-integer?", unOp isExactInteger),
    ("finite?", unOp isFinite),
    ("infinite?", unOp isInfinite'),
    ("nan?", unOp isNaN'),
    ("=", numericBoolOp (==)),
    ("<", numericBoolOp (<)),
    (">", numericBoolOp (>)),
    ("<=", numericBoolOp (<=)),
    (">=", numericBoolOp (>=)),
    ("+", add),
    ("-", sub),
    ("*", multiply),
    ("/", division),
    ("abs", unOp numAbs),
    ("floor/", floorInt),
    ("floor-quotient", floorQuot),
    ("floor-remainder", floorRem),
    ("truncate/", truncateInt),
    ("truncate-quotient", truncateQuot),
    ("truncate-remainder", truncateRem),
    ("quotient", truncateQuot),
    ("remainder", truncateRem),
    ("modulo", floorRem),
    ("gcd", numGcd),
    ("lcm", numLcm)
  ]

isNumber :: SchemeVal -> SchemeVal
isNumber (Number (Integer _)) = Boolean True
isNumber (Number (Real _)) = Boolean True
isNumber (Number (Rational _)) = Boolean True
isNumber (Number (Complex _)) = Boolean True
isNumber _ = Boolean False

isComplex :: SchemeVal -> SchemeVal
isComplex (Number (Integer _)) = Boolean True
isComplex (Number (Real _)) = Boolean True
isComplex (Number (Rational _)) = Boolean True
isComplex (Number (Complex _)) = Boolean True
isComplex _ = Boolean False

isReal :: SchemeVal -> SchemeVal
isReal (Number (Real _)) = Boolean True
isReal (Number (Rational _)) = Boolean True
isReal (Number (Integer _)) = Boolean True
isReal (Number (Complex x)) = if imagPart x == 0 then Boolean True else Boolean False
isReal _ = Boolean False

isRational :: SchemeVal -> SchemeVal
isRational (Number (Rational _)) = Boolean True
isRational (Number (Integer _)) = Boolean True
isRational (Number (Real x)) = Boolean (not $ isInfinite x)
isRational _ = Boolean False

isInteger :: SchemeVal -> SchemeVal
isInteger (Number (Integer _)) = Boolean True
isInteger (Number (Real x)) = Boolean (isDoubleInt x)
isInteger (Number (Rational x)) = Boolean (numerator x >= denominator x && numerator x `mod` denominator x == 0)
isInteger (Number (Complex x)) = Boolean (isDoubleInt (realPart x) && isDoubleInt (imagPart x))
isInteger _ = Boolean False

isExact :: SchemeVal -> SchemeVal
isExact (Number (Integer _)) = Boolean True
isExact (Number (Rational _)) = Boolean True
isExact _ = Boolean False

isInexact :: SchemeVal -> SchemeVal
isInexact (Number (Real _)) = Boolean True
isInexact _ = Boolean False

isExactInteger :: SchemeVal -> SchemeVal
isExactInteger x = case (isExact x, isInteger x) of
  (Boolean True, Boolean True) -> Boolean True
  _ -> Boolean False

isFinite :: SchemeVal -> SchemeVal
isFinite (Number (Real x)) = Boolean (not (isInfinite x || isNaN x))
isFinite (Number (Complex x)) = Boolean (not (isInfinite imag || isInfinite real || isNaN real || isNaN imag))
  where
    imag = imagPart x
    real = realPart x
isFinite _ = Boolean True

isInfinite' :: SchemeVal -> SchemeVal
isInfinite' (Number (Real x)) = Boolean (isInfinite x || isNaN x)
isInfinite' (Number (Complex x)) = Boolean (isInfinite imag || isInfinite real || isNaN real || isNaN imag)
  where
    imag = imagPart x
    real = realPart x
isInfinite' _ = Boolean False

isNaN' :: SchemeVal -> SchemeVal
isNaN' (Number (Real x)) = Boolean (isNaN x)
isNaN' (Number (Complex x)) = Boolean (isNaN (imagPart x) || isNaN (realPart x))
isNaN' _ = Boolean False

add :: [SchemeVal] -> SchemeResult SchemeVal
add [] = pure $ Number $ Integer 0
add xs = do
  nums <- mapM unwrapNumber xs
  return $ Number (sum nums)

multiply :: [SchemeVal] -> SchemeResult SchemeVal
multiply [] = pure $ Number $ Integer 1
multiply xs = do
  nums <- mapM unwrapNumber xs
  return $ Number (product nums)

sub :: [SchemeVal] -> SchemeResult SchemeVal
sub [] = throwError $ ArgumentLengthMismatch 1 []
sub [x] = do
  num <- unwrapNumber x
  return $ Number (negate num)
sub xs = do
  nums <- mapM unwrapNumber xs
  return $ Number (foldl1 (-) nums)

division :: [SchemeVal] -> SchemeResult SchemeVal
division [] = throwError $ ArgumentLengthMismatch 1 []
division [x] = do
  num <- unwrapNumber x
  return $ Number (recip num)
division xs = do
  nums <- mapM unwrapNumber xs
  return $ Number (foldl1 (/) nums)

numAbs :: SchemeVal -> SchemeVal
numAbs (Number n) = Number $ abs n
numAbs x = throw $ TypeMismatch "number" x

floorInt :: [SchemeVal] -> SchemeResult SchemeVal
floorInt [Number (Integer x), Number (Integer y)] = return $ List [Number $ Integer x', Number $ Integer y']
  where
    (x', y') = x `divMod` y
floorInt [x, Number (Integer _)] = throwError $ TypeMismatch "integer" x
floorInt [Number (Integer _), x] = throwError $ TypeMismatch "integer" x
floorInt _ = throwError $ ArgumentLengthMismatch 2 []

floorQuot :: [SchemeVal] -> SchemeResult SchemeVal
floorQuot [Number (Integer x), Number (Integer y)] = return $ Number $ Integer (x `div` y)
floorQuot [x, Number (Integer _)] = throwError $ TypeMismatch "integer" x
floorQuot [Number (Integer _), x] = throwError $ TypeMismatch "integer" x
floorQuot _ = throwError $ ArgumentLengthMismatch 2 []

floorRem :: [SchemeVal] -> SchemeResult SchemeVal
floorRem [Number (Integer x), Number (Integer y)] = return $ Number $ Integer (x `mod` y)
floorRem [x, Number (Integer _)] = throwError $ TypeMismatch "integer" x
floorRem [Number (Integer _), x] = throwError $ TypeMismatch "integer" x
floorRem _ = throwError $ ArgumentLengthMismatch 2 []

truncateInt :: [SchemeVal] -> SchemeResult SchemeVal
truncateInt [Number (Integer x), Number (Integer y)] = return $ List [Number $ Integer x', Number $ Integer y']
  where
    (x', y') = x `quotRem` y
truncateInt [x, Number (Integer _)] = throwError $ TypeMismatch "integer" x
truncateInt [Number (Integer _), x] = throwError $ TypeMismatch "integer" x
truncateInt _ = throwError $ ArgumentLengthMismatch 2 []

truncateQuot :: [SchemeVal] -> SchemeResult SchemeVal
truncateQuot [Number (Integer x), Number (Integer y)] = return $ Number $ Integer (x `quot` y)
truncateQuot [x, Number (Integer _)] = throwError $ TypeMismatch "integer" x
truncateQuot [Number (Integer _), x] = throwError $ TypeMismatch "integer" x
truncateQuot _ = throwError $ ArgumentLengthMismatch 2 []

truncateRem :: [SchemeVal] -> SchemeResult SchemeVal
truncateRem [Number (Integer x), Number (Integer y)] = return $ Number $ Integer (x `rem` y)
truncateRem [x, Number (Integer _)] = throwError $ TypeMismatch "integer" x
truncateRem [Number (Integer _), x] = throwError $ TypeMismatch "integer" x
truncateRem _ = throwError $ ArgumentLengthMismatch 2 []

numGcd :: [SchemeVal] -> SchemeResult SchemeVal
numGcd [] = return $ Number (Integer 0)
numGcd [x] = return x
numGcd xs = do
  ints <- mapM unwrapToInt xs
  return $ Number $ Integer (foldl1 gcd ints)

numLcm :: [SchemeVal] -> SchemeResult SchemeVal
numLcm [] = return $ Number (Integer 1)
numLcm [x] = return x
numLcm xs = do
  ints <- mapM unwrapToInt xs
  return $ Number $ Integer (foldl1 lcm ints)

numericBoolOp :: (Number -> Number -> Bool) -> [SchemeVal] -> Either SchemeError SchemeVal
numericBoolOp _ [] = throwError $ ArgumentLengthMismatch 1 []
numericBoolOp _ [x] = case isNaN' x of
  Boolean True -> return $ Boolean False
  _ -> return $ Boolean True
numericBoolOp op xs = do
  nums <- mapM unwrapNumber xs
  return $ Boolean $ and $ zipWith op nums $ drop 1 nums

isDoubleInt :: Double -> Bool
isDoubleInt d = (ceiling d :: Integer) == (floor d :: Integer)

unwrapToInt :: SchemeVal -> SchemeResult Integer
unwrapToInt (Number (Integer x)) = pure x
unwrapToInt (Number (Real x)) = pure $ floor x
unwrapToInt (Number (Rational x)) = pure $ floor x
unwrapToInt x = throwError $ TypeMismatch "integer" x

unwrapNumber :: SchemeVal -> SchemeResult Number
unwrapNumber (Number x) = pure x
unwrapNumber (List [x]) = unwrapNumber x
unwrapNumber x = throwError $ TypeMismatch "number" x

ext: hs


text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE ConstraintKinds       #-}
{-# LANGUAGE DeriveDataTypeable    #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
#if __GLASGOW_HASKELL__ >= 800
    {-# OPTIONS_GHC -fno-warn-redundant-constraints #-}
#endif
-- |
-- Module      : Graphics.Image.Interface.Repa.Generic
-- Copyright   : (c) Alexey Kuleshevich 2017
-- License     : BSD3
-- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
-- Stability   : experimental
-- Portability : non-portable
--
module Graphics.Image.Interface.Repa.Generic
  ( RImage
  , Strategy(..)
  , computeR
  , makeImageR
  , dimsR
  , scalarR
  , index00R
  , makeImageWindowedR
  , mapR
  , imapR
  , eqR
  , zipWithR
  , izipWithR
  , unsafeTraverseR
  , traverseR
  , unsafeTraverse2R
  , traverse2R
  , transposeR
  , backpermuteR
  , fromListsR
  , foldR
  , foldIxR
  , toVectorUnboxedR
  , fromVectorUnboxedR
  , toVectorStorableR
  , fromVectorStorableR
  , fromRepaArrayR
  , multR
  , ix2sh
  , sh2ix
  )where

import           Data.Array.Repa                     as R
import           Data.Array.Repa.Eval                as R (Elt (..), Target,
                                                           fromList,
                                                           suspendedComputeP)
import           Data.Array.Repa.Operators.Traversal as R
import           Data.Array.Repa.Repr.ForeignPtr
import           Data.Array.Repa.Repr.Partitioned    (P, Range (..))
import           Data.Array.Repa.Repr.Unboxed        (Unbox)
import           Data.Array.Repa.Repr.Undefined      (X)
import           Data.Complex                        as C
import           Data.Maybe                          (listToMaybe)
import qualified Data.Vector.Storable                as VS
import qualified Data.Vector.Unboxed                 as VU
import           Graphics.Image.ColorSpace.Binary    (Bit (..))
import           Graphics.Image.Interface            (ColorSpace (..), Pixel,
                                                      checkDims, toIx)
import           Prelude                             as P

data Strategy
  = Parallel
  | Sequential

data RImage r p = RTImage !(R.Array r R.DIM2 p)
                | RDImage !(R.Array R.D R.DIM2 p)

sh2ix :: DIM2 -> (Int, Int)
sh2ix (Z :. i :. j) = (i, j)
{-# INLINE sh2ix #-}

ix2sh :: (Int, Int) -> DIM2
ix2sh !(i, j) = Z :. i :. j
{-# INLINE ix2sh #-}

errorR :: String -> String -> a
errorR fName errMsg =
  error $ "Graphics.Image.Interface.Repa.Generic." P.++ fName P.++ ": " P.++ errMsg

checkDimsR :: String -> (Int, Int) -> DIM2
checkDimsR fName = ix2sh . checkDims ("Graphics.Image.Interface.Repa.Generic." P.++ fName)
{-# INLINE checkDimsR #-}

dimsR :: Source r p => RImage r p -> (Int, Int)
dimsR (RTImage arr) = sh2ix $ R.extent arr
dimsR (RDImage arr) = sh2ix $ R.extent arr
{-# INLINE dimsR #-}

makeImageR :: (Int, Int) -> ((Int, Int) -> px) -> RImage t px
makeImageR sz f =
  RDImage $ R.fromFunction (checkDimsR "makeImageR" sz) (f . sh2ix)
{-# INLINE makeImageR #-}

makeImageWindowedR
  :: (Int, Int)
  -> (Int, Int)
  -> (Int, Int)
  -> ((Int, Int) -> p)
  -> ((Int, Int) -> p)
  -> RImage r p
makeImageWindowedR !sz !wIx !wSz getWindowPx getBorderPx =
  RDImage $ R.delay $ makeWindowed sh wIx wSz wArr bArr
  where
    wArr = R.fromFunction sh (getWindowPx . sh2ix)
    bArr = R.fromFunction sh (getBorderPx . sh2ix)
    !sh = checkDimsR "makeImageWindowedR" sz
{-# INLINE makeImageWindowedR #-}

makeWindowed
  :: (Source r1 p, Source r2 p)
  => DIM2 -- ^ Extent of array.
  -> (Int, Int) -- ^ Window starting index
  -> (Int, Int) -- ^ Window size.
  -> Array r1 DIM2 p -- ^ Array for internal elements.
  -> Array r2 DIM2 p -- ^ Array for border elements.
  -> Array (P r1 (P r2 (P r2 (P r2 (P r2 X))))) DIM2 p
makeWindowed sh@(Z :. m :. n) !(it, jt) !(wm, wn) arrWindow arrBorder =
  let inInternal (Z :. i :. j) = i >= it && i < ib && j >= jt && j < jb
      {-# INLINE inInternal #-}
      inBorder = not . inInternal
      {-# INLINE inBorder #-}
      !(ib, jb) = (wm + it, wn + jt)
  in APart sh (Range (Z :. it :. jt) (Z :. wm       :. wn      ) inInternal) arrWindow $
     APart sh (Range (Z :. 0  :. 0 ) (Z :. it       :. n       ) inBorder  ) arrBorder $
     APart sh (Range (Z :. it :. 0 ) (Z :. wm       :. jt      ) inBorder  ) arrBorder $
     APart sh (Range (Z :. it :. jb) (Z :. wm       :. (n - jb)) inBorder  ) arrBorder $
     APart sh (Range (Z :. ib :. 0 ) (Z :. (m - ib) :. n       ) inBorder  ) arrBorder $
     AUndefined sh
{-# INLINE makeWindowed #-}

fromRepaArrayR :: (Source r1 p) => R.Array r1 DIM2 p -> RImage r p
fromRepaArrayR = RDImage . R.delay
{-# INLINE fromRepaArrayR #-}

scalarR :: p -> RImage t p
scalarR = RDImage . fromFunction (Z :. 1 :. 1) . const
{-# INLINE scalarR #-}

index00R :: Source r e => RImage r e -> e
index00R (RTImage arr) = R.index arr (Z :. 0 :. 0)
index00R (RDImage arr) = R.index arr (Z :. 0 :. 0)
{-# INLINE index00R #-}

mapR :: Source r1 a => (a -> p) -> RImage r1 a -> RImage r p
mapR f (RTImage arr) = RDImage (R.map f arr)
mapR f (RDImage arr) = RDImage (R.map f arr)
{-# INLINE mapR #-}

imapR :: Source r2 b =>
         ((Int, Int) -> b -> p) -> RImage r2 b -> RImage r p
imapR f (RTImage arr) = RDImage (imapArr (f . sh2ix) arr)
imapR f (RDImage arr) = RDImage (imapArr (f . sh2ix) arr)
{-# INLINE imapR #-}

imapArr :: R.Source r2 b =>
           (DIM2 -> b -> c) -> R.Array r2 DIM2 b -> R.Array R.D DIM2 c
imapArr f !arr = R.traverse arr id (\ getPx !sh -> f sh (getPx sh))
{-# INLINE imapArr #-}

getDelayed :: Source r e => RImage r e -> Array D DIM2 e
getDelayed (RTImage arr) = delay arr
getDelayed (RDImage arr) = arr
{-# INLINE getDelayed #-}

zipWithR
  :: (Source r2 t, Source r1 e) =>
     (t -> e -> c) -> RImage r2 t -> RImage r1 e -> RImage r c
zipWithR f !img1 !img2 = zipWithArr f (getDelayed img1) (getDelayed img2)
{-# INLINE zipWithR #-}

zipWithArr
  :: (Source r1 a, Source r2 b)
  => (a -> b -> c) -> Array r1 DIM2 a -> Array r2 DIM2 b -> RImage r c
zipWithArr f arr1 arr2 =
  RDImage $
  case (extent arr1, extent arr2) of
    (Z :. 1 :. 1, _) -> R.map (f (unsafeIndex arr1 (Z :. 0 :. 0))) arr2
    (_, Z :. 1 :. 1) -> R.map (`f` (unsafeIndex arr2 (Z :. 0 :. 0))) arr1
    _                -> R.zipWith f arr1 arr2
{-# INLINE zipWithArr #-}

izipWithR
  :: (Source r2 a, Source r1 b) =>
     ((Int, Int) -> a -> b -> c)
     -> RImage r2 a -> RImage r1 b -> RImage r c
izipWithR f !img1 !img2 = izipWithArr f (getDelayed img1) (getDelayed img2)
{-# INLINE izipWithR #-}

izipWithArr
  :: (Source r1 a, Source r2 b)
  => ((Int, Int) -> a -> b -> c)
  -> Array r1 DIM2 a
  -> Array r2 DIM2 b
  -> RImage r c
izipWithArr f arr1 arr2 =
  RDImage $
  case (extent arr1, extent arr2) of
    (Z :. 1 :. 1, _) ->
      imapArr (\ !sh -> f (sh2ix sh) (unsafeIndex arr1 (Z :. 0 :. 0))) arr2
    (_, Z :. 1 :. 1) ->
      imapArr
        (\ !sh !px1 -> f (sh2ix sh) px1 (unsafeIndex arr2 (Z :. 0 :. 0))) arr1
    (Z :. m1 :. n1, Z :. m2 :. n2) -> traverse2 arr1 arr2 getNewDims getNewPx
      where getNewDims _ _ = Z :. min m1 m2 :. min n1 n2
            {-# INLINE getNewDims #-}
            getNewPx getPx1 getPx2 !sh = f (sh2ix sh) (getPx1 sh) (getPx2 sh)
            {-# INLINE getNewPx #-}
{-# INLINE izipWithArr #-}

traverseR
  :: Source r1 c
  => RImage r1 c
  -> ((Int, Int) -> (Int, Int))
  -> (((Int, Int) -> c) -> (Int, Int) -> p)
  -> RImage r p
traverseR img getNewDims getNewPx =
  RDImage $
  R.traverse
    (getDelayed img)
    (checkDimsR "traverseR" . getNewDims . sh2ix)
    (\ getPx -> getNewPx (getPx . ix2sh) . sh2ix)
{-# INLINE traverseR #-}

unsafeTraverseR
  :: Source r1 c
  => RImage r1 c
  -> ((Int, Int) -> (Int, Int))
  -> (((Int, Int) -> c) -> (Int, Int) -> p)
  -> RImage r p
unsafeTraverseR img getNewDims getNewPx =
  RDImage $
  R.unsafeTraverse
    (getDelayed img)
    (checkDimsR "traverseR" . getNewDims . sh2ix)
    (\ getPx -> getNewPx (getPx . ix2sh) . sh2ix)
{-# INLINE unsafeTraverseR #-}

traverse2R
  :: (Source r1 a, Source r2 b)
  => RImage r1 a
  -> RImage r2 b
  -> ((Int, Int) -> (Int, Int) -> (Int, Int))
  -> (((Int, Int) -> a) -> ((Int, Int) -> b) -> (Int, Int) -> c)
  -> RImage r c
traverse2R img1 img2 getNewDims getNewPx =
  RDImage $
  R.traverse2
    (getDelayed img1)
    (getDelayed img2)
    (\ !sh1 !sh2 -> checkDimsR "traverse2R" (getNewDims (sh2ix sh1) (sh2ix sh2)))
    (\getPx1 getPx2 -> getNewPx (getPx1 . ix2sh) (getPx2 . ix2sh) . sh2ix)
{-# INLINE traverse2R #-}

unsafeTraverse2R
  :: (Source r1 a, Source r2 b)
  => RImage r1 a
  -> RImage r2 b
  -> ((Int, Int) -> (Int, Int) -> (Int, Int))
  -> (((Int, Int) -> a) -> ((Int, Int) -> b) -> (Int, Int) -> c)
  -> RImage r c
unsafeTraverse2R img1 img2 getNewDims getNewPx =
  RDImage $
  R.unsafeTraverse2
    (getDelayed img1)
    (getDelayed img2)
    (\ !sh1 !sh2 -> checkDimsR "traverse2R" (getNewDims (sh2ix sh1) (sh2ix sh2)))
    (\getPx1 getPx2 -> getNewPx (getPx1 . ix2sh) (getPx2 . ix2sh) . sh2ix)
{-# INLINE unsafeTraverse2R #-}

transposeR :: Source r1 p => RImage r1 p -> RImage r p
transposeR (RDImage arr) = RDImage (R.transpose arr)
transposeR (RTImage arr) = RDImage (R.transpose arr)
{-# INLINE transposeR #-}

backpermuteR
  :: R.Source r e
  => (Int, Int) -> ((Int, Int) -> (Int, Int)) -> RImage r e -> RImage r e
backpermuteR !newDims g !img =
  RDImage $ R.backpermute
    (ix2sh (checkDims "backpermuteR" newDims))
    (ix2sh . g . sh2ix)
    (getDelayed img)
{-# INLINE backpermuteR #-}

fromListsR :: Target r p => [[p]] -> RImage r p
fromListsR ls =
  if all (== n) (P.map length ls)
    then RTImage $ R.fromList (Z :. m :. n) . concat $ ls
    else error "fromListsRepa: Inner lists do not all have an equal length."
  where
    !(m, n) = checkDims "fromListsRepa" (length ls, maybe 0 length $ listToMaybe ls)
{-# INLINE fromListsR #-}

multR :: (Num p, Elt p, Unbox p, Source r2 p, Source r1 p, Source r p, Target r p,
           Target r2 p, Target r1 p) =>
         Strategy -> RImage r2 p -> RImage r1 p -> RImage r p
multR strategy (RTImage arr1) (RTImage arr2) =
  if n1 /= m2
    then errorR "multR" $
         "Inner dimensions of multiplied images must be the same, but received: " P.++
         show (m1, n1) P.++
         " X " P.++
         show (m2, n2)
    else computeR strategy $ RDImage $ fromFunction (Z :. m1 :. n2) getPx
  where
    (Z :. m1 :. n1) = R.extent arr1
    (Z :. m2 :. n2) = R.extent arr2
    getPx (Z :. i :. j) =
      R.sumAllS
        (R.slice arr1 (R.Any :. (i :: Int) :. R.All) R.*^
         R.slice arr2 (R.Any :. (j :: Int)))
    {-# INLINE getPx #-}
multR strategy img1 img2 =
  multR strategy (computeR strategy img1) (computeR strategy img2)
{-# INLINE multR #-}

computeR :: (Target r p, Source r p) => Strategy -> RImage r p -> RImage r p
computeR Parallel (RDImage arr) = arrManifest `deepSeqArray` RTImage arrManifest
  where arrManifest = suspendedComputeP arr
computeR Sequential (RDImage arr) = RTImage $ computeS arr
computeR _ img = img
{-# INLINE computeR #-}

foldR :: (Source r p, Target r p, Elt p, Unbox p)
       => Strategy -> (p -> p -> p) -> p -> RImage r p -> p
foldR Parallel f !px0 !img =
  let RTImage arr = computeR Parallel img
  in case foldAllP f px0 arr of
       Just e  -> e
       Nothing -> errorR "foldPR" "impossible happened."
foldR Sequential f !px0 img =
  let RTImage arr = computeR Sequential img
  in foldAllS f px0 arr
{-# INLINE foldR #-}

addIxArr
  :: Source r2 b =>
     Array r2 DIM2 b -> Array R.D DIM2 (Int, b)
addIxArr !arr = R.zipWith (,) arrIx arr
  where
    !sh = extent arr
    !arrIx = R.fromFunction (R.extent arr) (toIndex sh)
{-# INLINE addIxArr #-}

foldIxR :: (Source r b, Elt b, Unbox b) =>
           Strategy -> (b -> (Int, Int) -> b -> b) -> b -> RImage r b -> b
foldIxR strategy f !acc' !img =
  case strategy of
    Parallel ->
      case foldAllP accumWithIx (-1, acc') arr of
        Just (_, acc) -> acc
        Nothing       -> error $ "foldIxPR: impossible happened."
    Sequential -> snd $ R.foldAllS accumWithIx (-1, acc') arr
  where
    (Z :. _ :. n) = extent arr
    !arr = addIxArr (getDelayed img)
    accumWithIx !(-1, acc) !(k, px) = (-1, f acc (toIx n k) px)
    accumWithIx !(k, px) (-1, acc) = (-1, f acc (toIx n k) px)
    accumWithIx (acc1Ix, _) (acc2Ix, _) =
      errorR "accumWithx" $
      "Impossible happened. Received: " P.++ show acc1Ix P.++ " " P.++
      show acc2Ix
    {-# INLINE accumWithIx #-}
{-# INLINE foldIxR #-}

eqR :: (Eq a, Source r1 a, Source r a) =>
        Strategy -> RImage r1 a -> RImage r a -> Bool
eqR strategy !img1 !img2 =
  case strategy of
    Parallel ->
      let eqArr = R.zipWith (==) (getDelayed img1) (getDelayed img2)
      in dimsR img1 == dimsR img2 && mCompute "eqR" (foldAllP (&&) True eqArr)
    Sequential -> R.equalsS (getDelayed img1) (getDelayed img2)
{-# INLINE eqR #-}

mCompute :: String -> Maybe t -> t
mCompute _ (Just res)  = res
mCompute fName Nothing = errorR fName "impossible happened"
{-# INLINE mCompute #-}

fromVectorUnboxedR :: Unbox e => (Int, Int) -> VU.Vector e -> RImage U e
fromVectorUnboxedR !sz = RTImage . fromUnboxed (ix2sh sz)
{-# INLINE fromVectorUnboxedR #-}

toVectorR :: (Target r2 e, Source r2 e) =>
             (R.Array r2 DIM2 e -> t) -> Strategy -> RImage r2 e -> t
toVectorR toVector Parallel (RDImage arr) =
 toVector $ mCompute "toVectorR" (computeP arr)
toVectorR toVector Sequential (RDImage arr) = toVector $ computeS arr
toVectorR toVector _ (RTImage arr) = toVector arr
{-# INLINE toVectorR #-}

toVectorUnboxedR :: Unbox e => Strategy -> RImage U e -> VU.Vector e
toVectorUnboxedR = toVectorR toUnboxed
{-# INLINE toVectorUnboxedR #-}

toVectorStorableR :: VS.Storable e => Strategy -> RImage F e -> VS.Vector e
toVectorStorableR = toVectorR toStorableR
{-# INLINE toVectorStorableR #-}

toStorableR :: VS.Storable a => Array F DIM2 a -> VS.Vector a
toStorableR arr = VS.unsafeFromForeignPtr0 (toForeignPtr arr) (m * n)
  where
    (Z :. m :. n) = R.extent arr
{-# INLINE toStorableR #-}

fromVectorStorableR
  :: VS.Storable px
  => (Int, Int) -> VS.Vector px -> RImage F px
fromVectorStorableR !(m, n) !v
  | sz == sz' = RTImage $ fromForeignPtr (ix2sh (m, n)) fp
  | otherwise =
    errorR "fromVectorStorableR" $
    "(impossible) Vector size mismatch: " P.++ show sz P.++ " vs " P.++ show sz'
  where
    !(fp, sz) = VS.unsafeToForeignPtr0 v
    !sz' = m * n
{-# INLINE fromVectorStorableR #-}

instance (Num e, R.Elt e) => Elt (C.Complex e) where
  touch (r :+ i) = R.touch r >> R.touch i
  {-# INLINE touch #-}

  zero     = 0 :+ 0
  {-# INLINE zero #-}

  one      = 1 :+ 0
  {-# INLINE one #-}

instance Elt Bit where
  touch (Bit w) = R.touch w
  {-# INLINE touch #-}

  zero     = 0
  {-# INLINE zero #-}

  one      = 1
  {-# INLINE one #-}

instance (ColorSpace cs e, Elt e, Num (Pixel cs e)) => Elt (Pixel cs e) where
  touch = foldlPx (\ a !e -> a >> touch e) (return ())
  {-# INLINE touch #-}

  zero     = 0
  {-# INLINE zero #-}

  one      = 1
  {-# INLINE one #-}

max_stars_repo_name: khilanravani/hip

max_issues_repo_name: khilanravani/hip

max_forks_repo_name: khilanravani/hip


text:
{-# LANGUAGE FlexibleContexts #-}

import           AI.Layer
import           AI.Neuron

import           AI.Network
import           AI.Network.FeedForwardNetwork

import           AI.Trainer
import           AI.Trainer.BackpropTrainer

--import Network.Visualizations
import           Numeric.LinearAlgebra
import           System.IO
import           System.Random

main :: IO ()
main = do

  g <- newStdGen
  let l = LayerDefinition sigmoidNeuron 2 connectFully randomizeFully
  let l' = LayerDefinition sigmoidNeuron 2 connectFully randomizeFully
  let l'' = LayerDefinition sigmoidNeuron 1 connectFully randomizeFully

  let n = createNetwork normals g [l, l', l'']

  let t = BackpropTrainer (3 :: Double) quadraticCost quadraticCost'
  let dat = [(fromList [0, 1], fromList [1]), (fromList [1, 1], fromList [0]), (fromList [1, 0], fromList [1]), (fromList [0, 0], fromList [0])]

  let n' = trainNTimes g n t online dat 1000

  putStrLn "==> XOR predictions: "
  print $ predict (fromList [0, 0]) n'
  print $ predict (fromList [1, 0]) n'
  print $ predict (fromList [0, 1]) n'
  print $ predict (fromList [1, 1]) n'

  saveFeedForwardNetwork "xor.ann" n'

  putStrLn "==> Network saved and reloaded: "
  n'' <- loadFeedForwardNetwork "xor.ann" [l, l', l'']

  print $ predict (fromList [0, 0]) n''
  print $ predict (fromList [1, 0]) n''
  print $ predict (fromList [0, 1]) n''
  print $ predict (fromList [1, 1]) n''

  --networkHistogram "weights.png" weightList n''
  --networkHistogram "biases.png" biasList n''

max_stars_repo_name: jbarrow/LambdaNet

max_issues_repo_name: jbarrow/LambdaNet

max_forks_repo_name: jbarrow/LambdaNet


text:
{-# LANGUAGE NoMonomorphismRestriction #-}

module Main where

import LinAlg (testCols, testRows)

import Numeric.LinearAlgebra.Data
import Numeric.LinearAlgebra.HMatrix

import Data.Vector.Storable

import Prelude hiding (sum)

testMult n m = sum $ mat #> vec
  where
    mat = matrix n $ fmap fromIntegral [1..(n*m)]
    vec = fromList $ fmap fromIntegral [1..n]

main = print $ testMult testCols testRows


text:
module Problem215 where
{-
  Task description:
  Consider the problem of building a wall out of 2×1 and 3×1 bricks (horizontal×vertical dimensions) such that,
  for extra strength, the gaps between horizontally-adjacent bricks never line up in consecutive layers, i.e. never form a "running crack".

  For example, the following 9×3 wall is not acceptable due to the running crack shown in red:

  There are eight ways of forming a crack-free 9×3 wall, written W(9,3) = 8.

  Calculate W(32,10).
-}

import Control.Monad
import Data.Function (on)
import Data.Graph (Graph, Vertex)
import Data.Maybe (fromJust)
import Numeric.LinearAlgebra (Matrix, (<>))
import qualified GHC.Arr as Arr
import qualified Data.Graph as Graph
import qualified Data.HashSet as HashSet
import qualified Data.List as List
import qualified Data.Map as Map
import qualified Data.Tree as Tree
import qualified Numeric.LinearAlgebra as LA

{-
A row may only ever depend on the preceding row.
We may have cycles of rows that are allowed to follow each other.
An interesting question will be wether collapsing such cycles will be easier than simply 'exploring them'.
In both ways memoization may come in handy here.
-}

type Brick = Int
type Line = [Brick]
type Width = Int
type Height = Int
type Wall = [Line]

bricks :: [Brick]
bricks = [2, 3]

buildLines :: Width -> [Line]
buildLines 0 = [[]]
buildLines w = do
  brick <- bricks
  guard $ brick <= w
  lineRest <- buildLines $ w - brick
  return $ brick:lineRest

matchingLines :: Line -> Line -> Bool
matchingLines l1 l2 = not $ (matchingCracks `on` crackPositions) l1 l2
  where
    crackPositions :: Line -> [Int]
    crackPositions = fmap sum . init . tail . List.inits

    matchingCracks :: [Int] -> [Int] -> Bool
    matchingCracks [] _ = False
    matchingCracks _ [] = False
    matchingCracks xss@(x:xs) yss@(y:ys)
      | x == y = True
      | x < y = matchingCracks xs yss
      | x > y = matchingCracks xss ys

{-
With the matchingLines predicate I can generate a graph out of the possible lines.
We're searching for all routes of a given length that traverse this graph
and that start in any of it's nodes.
-}

mkLineEdges :: [Line] -> [(Line, Line)]
mkLineEdges [] = []
mkLineEdges (l:ls) = findMatchingLines l ls `mappend` mkLineEdges ls
  where
    findMatchingLines :: Line -> [Line] -> [(Line, Line)]
    findMatchingLines l lines = concat $ do
      l' <- lines
      guard $ matchingLines l l'
      return [(l,l'), (l',l)]

mkLineToVertex :: [Line] -> (Line -> Vertex)
mkLineToVertex lines = let lvMap = Map.fromList $ zip lines [0..]
                       in fromJust . (`Map.lookup` lvMap)

graphForWidth :: Width -> Graph
graphForWidth w = let lines = buildLines w
                      lEdges = mkLineEdges lines
                      bounds = (0, length lines - 1)
                      toVertex = mkLineToVertex lines
                      edges = [(toVertex a, toVertex b)|(a,b) <- lEdges]
                  in Graph.buildG bounds edges

graphToDenseMatrix :: Graph -> Matrix Double
graphToDenseMatrix = LA.toDense . graphToAssocMatrix
  where
    graphToAssocMatrix :: Graph -> LA.AssocMatrix
    graphToAssocMatrix g = do
      (v, targets) <- Arr.assocs g
      [((v, t), 1)| t <- targets]

dissect :: Graph -> [Graph]
dissect g = do
  tree <- Graph.components g
  let vertices = Tree.flatten tree
      (renumber, bounds) = mkRenumber vertices
  return . Graph.buildG bounds $ do
    v <- vertices
    w <-  Arr.unsafeAt g v
    return (renumber v, renumber w)
  where
    mkRenumber ls = let mapping = Map.fromList $ zip ls [0..]
                        λ = fromJust . (`Map.lookup` mapping)
                    in (λ, (0, Map.size mapping))

mPow :: Matrix Double -> Int -> Matrix Double
mPow m 0 = LA.ident $ LA.rows m
mPow m 1 = m
mPow m n
  | even n = mPow (m <> m) (n `div` 2)
  | otherwise = m <> mPow m (n - 1)

countWalls :: Width -> Height -> Int
countWalls w h
  | h < 1 = 0
  | h == 1 = length $ buildLines w
-- Using dissect reduced computation time from ~140s to ~9s.
  | otherwise = let g = graphForWidth w
                    ms = graphToDenseMatrix <$> dissect g
                    ms' = fmap (mPow `flip` (h - 1)) ms
                    calcSum = sum . concat . LA.toLists
                in round . sum $ fmap calcSum ms'

example = countWalls 9 3

-- 806844323190414
solution = countWalls 32 10

main = print solution

max_stars_repo_name: runjak/projectEuler

max_issues_repo_name: runjak/projectEuler

max_forks_repo_name: runjak/projectEuler


text:
{-# LANGUAGE CPP                    #-}
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE ScopedTypeVariables    #-}
{-# LANGUAGE ViewPatterns           #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

-----------------------------------------------------------------------------
{- |
Module      :  Internal.Util
Copyright   :  (c) Alberto Ruiz 2013
License     :  BSD3
Maintainer  :  Alberto Ruiz
Stability   :  provisional

-}
-----------------------------------------------------------------------------

module Internal.Util(

    -- * Convenience functions
    vector, matrix,
    disp,
    formatSparse,
    approxInt,
    dispDots,
    dispBlanks,
    formatShort,
    dispShort,
    zeros, ones,
    diagl,
    row,
    col,
    (&), (¦), (|||), (——), (===),
    (?), (¿),
    Indexable(..), size,
    Numeric,
    rand, randn,
    cross,
    norm,
    ℕ,ℤ,ℝ,ℂ,iC,
    Normed(..), norm_Frob, norm_nuclear,
    magnit,
    normalize,
    mt,
    (~!~),
    pairwiseD2,
    rowOuters,
    null1,
    null1sym,
    -- * Convolution
    -- ** 1D
    corr, conv, corrMin,
    -- ** 2D
    corr2, conv2, separable,
    block2x2,block3x3,view1,unView1,foldMatrix,
    gaussElim_1, gaussElim_2, gaussElim,
    luST, luSolve', luSolve'', luPacked', luPacked'',
    invershur
) where

import           Control.Arrow        ((&&&), (***))
import           Control.Monad        (forM_, when)
import           Data.Complex
import           Data.Function        (on)
import           Data.List            (foldl', intercalate, sortBy)
import           Data.List.Split      (splitOn)
import           Internal.Algorithms  hiding (Normed, linearSolve', luPacked', luSolve')
import           Internal.Container
import           Internal.Convolution
import           Internal.Element
import           Internal.IO
import           Internal.Matrix      hiding (size)
import           Internal.Numeric
import           Internal.Random
import           Internal.ST
import           Internal.Vector
import           Internal.Vectorized
import           Numeric.Matrix       ()
import           Numeric.Vector       ()
import           Text.Printf
#if MIN_VERSION_base(4,11,0)
import           Prelude              hiding ((<>))
#endif

type ℝ = Float
type ℕ = Int
type ℤ = Int
type ℂ = Complex Float

-- | imaginary unit
iC :: C
iC = 0:+1

{- | Create a real vector.

>>> vector [1..5]
[1.0,2.0,3.0,4.0,5.0]
it :: Vector R

-}
vector :: [R] -> Vector R
vector = fromList

{- | Create a real matrix.

>>> matrix 5 [1..15]
(3><5)
 [  1.0,  2.0,  3.0,  4.0,  5.0
 ,  6.0,  7.0,  8.0,  9.0, 10.0
 , 11.0, 12.0, 13.0, 14.0, 15.0 ]

-}
matrix
  :: Int -- ^ number of columns
  -> [R] -- ^ elements in row order
  -> Matrix R
matrix c = reshape c . fromList

{- | print a real matrix with given number of digits after the decimal point

>>> disp 5 $ ident 2 / 3
2x2
0.33333  0.00000
0.00000  0.33333

-}
disp :: Int -> Matrix Float -> IO ()

disp n = putStr . dispf n

{- | create a real diagonal matrix from a list

>>> diagl [1,2,3]
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 2.0, 0.0
 , 0.0, 0.0, 3.0 ]

-}
diagl :: [Float] -> Matrix Float
diagl = diag . fromList

-- | a real matrix of zeros
zeros :: Int -- ^ rows
      -> Int -- ^ columns
      -> Matrix Float
zeros r c = konst 0 (r,c)

-- | a real matrix of ones
ones :: Int -- ^ rows
     -> Int -- ^ columns
     -> Matrix Float
ones r c = konst 1 (r,c)

-- | concatenation of real vectors
infixl 3 &
(&) :: Vector Float -> Vector Float -> Vector Float
a & b = vjoin [a,b]

{- | horizontal concatenation

>>> ident 3 ||| konst 7 (3,4)
(3><7)
 [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0
 , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0
 , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ]

-}
infixl 3 |||
(|||) :: Element t => Matrix t -> Matrix t -> Matrix t
a ||| b = fromBlocks [[a,b]]

-- | a synonym for ('|||') (unicode 0x00a6, broken bar)
infixl 3 ¦
(¦) :: Matrix Float -> Matrix Float -> Matrix Float
(¦) = (|||)

-- | vertical concatenation
--
(===) :: Element t => Matrix t -> Matrix t -> Matrix t
infixl 2 ===
a === b = fromBlocks [[a],[b]]

-- | a synonym for ('===') (unicode 0x2014, em dash)
(——) :: Matrix Float -> Matrix Float -> Matrix Float
infixl 2 ——
(——) = (===)

-- | create a single row real matrix from a list
--
-- >>> row [2,3,1,8]
-- (1><4)
--  [ 2.0, 3.0, 1.0, 8.0 ]
--
row :: [Float] -> Matrix Float
row = asRow . fromList

-- | create a single column real matrix from a list
--
-- >>> col [7,-2,4]
-- (3><1)
--  [  7.0
--  , -2.0
--  ,  4.0 ]
--
col :: [Float] -> Matrix Float
col = asColumn . fromList

{- | extract rows

>>> (20><4) [1..] ? [2,1,1]
(3><4)
 [ 9.0, 10.0, 11.0, 12.0
 , 5.0,  6.0,  7.0,  8.0
 , 5.0,  6.0,  7.0,  8.0 ]

-}
infixl 9 ?
(?) :: Element t => Matrix t -> [Int] -> Matrix t
(?) = flip extractRows

{- | extract columns

(unicode 0x00bf, inverted question mark, Alt-Gr ?)

>>> (3><4) [1..] ¿ [3,0]
(3><2)
 [  4.0, 1.0
 ,  8.0, 5.0
 , 12.0, 9.0 ]

-}
infixl 9 ¿
(¿) :: Element t => Matrix t -> [Int] -> Matrix t
(¿)= flip extractColumns

cross :: Product t => Vector t -> Vector t -> Vector t
-- ^ cross product (for three-element vectors)
cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3]
          | otherwise = error $ "the cross product requires 3-element vectors (sizes given: "
                                ++show (dim x)++" and "++show (dim y)++")"
  where
    [x1,x2,x3] = toList x
    [y1,y2,y3] = toList y
    z1 = x2*y3-x3*y2
    z2 = x3*y1-x1*y3
    z3 = x1*y2-x2*y1

{-# SPECIALIZE cross :: Vector Float -> Vector Float -> Vector Float #-}
{-# SPECIALIZE cross :: Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float) #-}

norm :: Vector Float -> Float
-- ^ 2-norm of real vector
norm = pnorm PNorm2

-- | p-norm for vectors, operator norm for matrices
class Normed a
  where
    norm_0   :: a -> R
    norm_1   :: a -> R
    norm_2   :: a -> R
    norm_Inf :: a -> R

instance Normed (Vector R)
  where
    norm_0 v = sumElements (step (abs v - scalar (eps*normInf v)))
    norm_1 = pnorm PNorm1
    norm_2 = pnorm PNorm2
    norm_Inf = pnorm Infinity

instance Normed (Vector C)
  where
    norm_0 v = sumElements (step (fst (fromComplex (abs v)) - scalar (eps*normInf v)))
    norm_1 = pnorm PNorm1
    norm_2 = pnorm PNorm2
    norm_Inf = pnorm Infinity

instance Normed (Matrix R)
  where
    norm_0 = norm_0 . flatten
    norm_1 = pnorm PNorm1
    norm_2 = pnorm PNorm2
    norm_Inf = pnorm Infinity

instance Normed (Matrix C)
  where
    norm_0 = norm_0 . flatten
    norm_1 = pnorm PNorm1
    norm_2 = pnorm PNorm2
    norm_Inf = pnorm Infinity

instance Normed (Vector I)
  where
    norm_0 = fromIntegral . sumElements . step . abs
    norm_1 = fromIntegral . norm1
    norm_2 v = sqrt . fromIntegral $ dot v v
    norm_Inf = fromIntegral . normInf

instance Normed (Vector Z)
  where
    norm_0 = fromIntegral . sumElements . step . abs
    norm_1 = fromIntegral . norm1
    norm_2 v = sqrt . fromIntegral $ dot v v
    norm_Inf = fromIntegral . normInf

-- instance Normed (Vector Float)
--   where
--     norm_0 = norm_0 . double
--     norm_1 = norm_1 . double
--     norm_2 = norm_2 . double
--     norm_Inf = norm_Inf . double

-- instance Normed (Vector (Complex Float))
--   where
--     norm_0 = norm_0 . double
--     norm_1 = norm_1 . double
--     norm_2 = norm_2 . double
--     norm_Inf = norm_Inf . double

-- | Frobenius norm (Schatten p-norm with p=2)
norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R
norm_Frob = norm_2 . flatten

-- | Sum of singular values (Schatten p-norm with p=1)
norm_nuclear :: Field t => Matrix t -> R
norm_nuclear = sumElements . singularValues

{- | Check if the absolute value or complex magnitude is greater than a given threshold

>>> magnit 1E-6 (1E-12 :: R)
False
>>> magnit 1E-6 (3+iC :: C)
True
>>> magnit 0 (3 :: I ./. 5)
True

-}
magnit :: (Element t, Normed (Vector t)) => R -> t -> Bool
magnit e x = norm_1 (fromList [x]) > e

-- | Obtains a vector in the same direction with 2-norm=1
normalize :: (Normed (Vector t), Num (Vector t), Field t) => Vector t -> Vector t
normalize v = v / real (scalar (norm_2 v))

-- | trans . inv
mt :: Matrix Float -> Matrix Float
mt = trans . inv

--------------------------------------------------------------------------------
{- |

>>> size $ vector [1..10]
10
>>> size $ (2><5)[1..10::Float]
(2,5)

-}
size :: Container c t => c t -> IndexOf c
size = size'

{- | Alternative indexing function.

>>> vector [1..10] ! 3
4.0

On a matrix it gets the k-th row as a vector:

>>> matrix 5 [1..15] ! 1
[6.0,7.0,8.0,9.0,10.0]
it :: Vector Float

>>> matrix 5 [1..15] ! 1 ! 3
9.0

-}
class Indexable c t | c -> t , t -> c
  where
    infixl 9 !
    (!) :: c -> Int -> t

instance Indexable (Vector Float) Float
  where
    (!) = (@>)

-- instance Indexable (Vector Float) Float
--   where
--     (!) = (@>)

instance Indexable (Vector I) I
  where
    (!) = (@>)

instance Indexable (Vector Z) Z
  where
    (!) = (@>)

-- instance Indexable (Vector (Complex Float)) (Complex Float)
--   where
--     (!) = (@>)

instance Indexable (Vector (Complex Float)) (Complex Float)
  where
    (!) = (@>)

instance Element t => Indexable (Matrix t) (Vector t)
  where
    (!) m j = subVector (j*c) c (flatten m)
      where
        c = cols m

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

-- | Matrix of pairwise squared distances of row vectors
-- (using the matrix product trick in blog.smola.org)
pairwiseD2 :: Matrix Float -> Matrix Float -> Matrix Float
pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y
               | otherwise = error $ "pairwiseD2 with different number of columns: "
                                   ++ show (size x) ++ ", " ++ show (size y)
  where
    ox = one (rows x)
    oy = one (rows y)
    oc = one (cols x)
    one k = konst 1 k
    x2 = x * x <> oc
    y2 = y * y <> oc
    ok = cols x == cols y

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

{- | outer products of rows

>>> a
(3><2)
 [   1.0,   2.0
 ,  10.0,  20.0
 , 100.0, 200.0 ]
>>> b
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]

>>> rowOuters a (b ||| 1)
(3><8)
 [   1.0,   2.0,   3.0,   1.0,    2.0,    4.0,    6.0,   2.0
 ,  40.0,  50.0,  60.0,  10.0,   80.0,  100.0,  120.0,  20.0
 , 700.0, 800.0, 900.0, 100.0, 1400.0, 1600.0, 1800.0, 200.0 ]

-}
rowOuters :: Matrix Float -> Matrix Float -> Matrix Float
rowOuters a b = a' * b'
  where
    a' = kronecker a (ones 1 (cols b))
    b' = kronecker (ones 1 (cols a)) b

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

-- | solution of overconstrained homogeneous linear system
null1 :: Matrix R -> Vector R
null1 = last . toColumns . snd . rightSV

-- | solution of overconstrained homogeneous symmetric linear system
null1sym :: Herm R -> Vector R
null1sym = last . toColumns . snd . eigSH

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

infixl 0 ~!~
c ~!~ msg = when c (error msg)

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

formatSparse :: String -> String -> String -> Int -> Matrix Float -> String

formatSparse zeroI _zeroF sep _ (approxInt -> Just m) = format sep f m
  where
    f 0 = zeroI
    f x = printf "%.0f" x

formatSparse zeroI zeroF sep n m = format sep f m
  where
    f x | abs (x::Float) < 2*peps = zeroI++zeroF
        | abs (fromIntegral (round x::Int) - x) / abs x < 2*peps
            = printf ("%.0f."++replicate n ' ') x
        | otherwise = printf ("%."++show n++"f") x

approxInt m
    | norm_Inf (v - vi) < 2*peps * norm_Inf v = Just (reshape (cols m) vi)
    | otherwise = Nothing
  where
    v = flatten m
    vi = roundVector v

dispDots n = putStr . formatSparse "." (replicate n ' ') "  " n

dispBlanks n = putStr . formatSparse "" "" "  " n

formatShort sep fmt maxr maxc m = auxm4
  where
    (rm,cm) = size m
    (r1,r2,r3)
        | rm <= maxr = (rm,0,0)
        | otherwise  = (maxr-3,rm-maxr+1,2)
    (c1,c2,c3)
        | cm <= maxc = (cm,0,0)
        | otherwise  = (maxc-3,cm-maxc+1,2)
    [ [a,_,b]
     ,[_,_,_]
     ,[c,_,d]] = toBlocks [r1,r2,r3]
                          [c1,c2,c3] m
    auxm = fromBlocks [[a,b],[c,d]]
    auxm2
        | cm > maxc = format "|" fmt auxm
        | otherwise = format sep fmt auxm
    auxm3
        | cm > maxc = map (f . splitOn "|") (lines auxm2)
        | otherwise = (lines auxm2)
    f items = intercalate sep (take (maxc-3) items) ++ "  .. " ++
              intercalate sep (drop (maxc-3) items)
    auxm4
        | rm > maxr = unlines (take (maxr-3) auxm3 ++ vsep : drop (maxr-3) auxm3)
        | otherwise = unlines auxm3
    vsep = map g (head auxm3)
    g '.' = ':'
    g _   = ' '

dispShort :: Int -> Int -> Int -> Matrix Float -> IO ()
dispShort maxr maxc dec m =
    printf "%dx%d\n%s" (rows m) (cols m) (formatShort "  " fmt maxr maxc m)
  where
    fmt = printf ("%."++show dec ++"f")

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

-- matrix views

block2x2 r c m = [[m11,m12],[m21,m22]]
  where
    m11 = m ?? (Take r, Take c)
    m12 = m ?? (Take r, Drop c)
    m21 = m ?? (Drop r, Take c)
    m22 = m ?? (Drop r, Drop c)

block3x3 r nr c nc m = [[m ?? (er !! i, ec !! j) | j <- [0..2] ] | i <- [0..2] ]
  where
    er = [ Range 0 1 (r-1), Range r 1 (r+nr-1), Drop (nr+r) ]
    ec = [ Range 0 1 (c-1), Range c 1 (c+nc-1), Drop (nc+c) ]

view1 :: Numeric t => Matrix t -> Maybe (View1 t)
view1 m
    | rows m > 0 && cols m > 0 = Just (e, flatten m12, flatten m21 , m22)
    | otherwise = Nothing
  where
    [[m11,m12],[m21,m22]] = block2x2 1 1 m
    e = m11 `atIndex` (0, 0)

unView1 :: Numeric t => View1 t -> Matrix t
unView1 (e,r,c,m) = fromBlocks [[scalar e, asRow r],[asColumn c, m]]

type View1 t = (t, Vector t, Vector t, Matrix t)

foldMatrix :: Numeric t => (Matrix t -> Matrix t) -> (View1 t -> View1 t) -> (Matrix t -> Matrix t)
foldMatrix g f ( (f <$>) . view1 . g -> Just (e,r,c,m)) = unView1 (e, r, c, foldMatrix g f m)
foldMatrix _ _ m = m

swapMax k m
    | rows m > 0 && j>0 = (j, m ?? (Pos (idxs swapped), All))
    | otherwise  = (0,m)
  where
    j = maxIndex $ abs (tr m ! k)
    swapped = j:[1..j-1] ++ 0:[j+1..rows m-1]

down g a = foldMatrix g f a
  where
    f (e,r,c,m)
        | e /= 0    = (1, r', 0, m - outer c r')
        | otherwise = error "singular!"
      where
        r' = r / scalar e

-- | generic reference implementation of gaussian elimination
--
-- @a <> gaussElim a b = b@
--
gaussElim_2
  :: (Eq t, Fractional t, Num (Vector t), Numeric t)
  => Matrix t -> Matrix t -> Matrix t

gaussElim_2 a b = flipudrl r
  where
    flipudrl = flipud . fliprl
    splitColsAt n = (takeColumns n &&& dropColumns n)
    go f x y = splitColsAt (cols a) (down f $ x ||| y)
    (a1,b1) = go (snd . swapMax 0) a b
    ( _, r) = go id (flipudrl $ a1) (flipudrl $ b1)

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

gaussElim_1
  :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
  => Matrix t -> Matrix t -> Matrix t

gaussElim_1 x y = dropColumns (rows x) (flipud $ fromRows s2)
  where
    rs = toRows $ x ||| y
    s1 = fromRows $ pivotDown (rows x) 0 rs      -- interesting
    s2 = pivotUp (rows x-1) (toRows $ flipud s1)

pivotDown
  :: forall t . (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
  => Int -> Int -> [Vector t] -> [Vector t]
pivotDown t n xs
    | t == n    = []
    | otherwise = y : pivotDown t (n+1) ys
  where
    y:ys = redu (pivot n xs)

    pivot k = (const k &&& id)
            . sortBy (flip compare `on` (abs. (!k)))

    redu :: (Int, [Vector t]) -> [Vector t]
    redu (k,x:zs)
        | p == 0 = error "gauss: singular!"  -- FIXME
        | otherwise = u : map f zs
      where
        p = x!k
        u = scale (recip (x!k)) x
        f z = z - scale (z!k) u
    redu (_,[]) = []

pivotUp
  :: forall t . (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
  => Int -> [Vector t] -> [Vector t]
pivotUp n xs
    | n == -1 = []
    | otherwise = y : pivotUp (n-1) ys
  where
    y:ys = redu' (n,xs)

    redu' :: (Int, [Vector t]) -> [Vector t]
    redu' (k,x:zs) = u : map f zs
      where
        u = x
        f z = z - scale (z!k) u
    redu' (_,[]) = []

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

gaussElim a b = dropColumns (rows a) $ fst $ mutable gaussST (a ||| b)

gaussST (r,_) x = do
    let n = r-1
        axpy m a i j = rowOper (AXPY a i j AllCols)     m
        swap m i j   = rowOper (SWAP i j AllCols)       m
        scal m a i   = rowOper (SCAL a (Row i) AllCols) m
    forM_ [0..n] $ \i -> do
        c <- maxIndex . abs . flatten <$> extractMatrix x (FromRow i) (Col i)
        swap x i (i+c)
        a <- readMatrix x i i
        when (a == 0) $ error "singular!"
        scal x (recip a) i
        forM_ [i+1..n] $ \j -> do
            b <- readMatrix x j i
            axpy x (-b) i j
    forM_ [n,n-1..1] $ \i -> do
        forM_ [i-1,i-2..0] $ \j -> do
            b <- readMatrix x j i
            axpy x (-b) i j

luST ok (r,_) x = do
    let axpy m a i j = rowOper (AXPY a i j (FromCol (i+1))) m
        swap m i j   = rowOper (SWAP i j AllCols)           m
    p <- newUndefinedVector r
    forM_ [0..r-1] $ \i -> do
        k <- maxIndex . abs . flatten <$> extractMatrix x (FromRow i) (Col i)
        writeVector p i (k+i)
        swap x i (i+k)
        a <- readMatrix x i i
        when (ok a) $ do
            forM_ [i+1..r-1] $ \j -> do
                b <- (/a) <$> readMatrix x j i
                axpy x (-b) i j
                writeMatrix x j i b
    v <- unsafeFreezeVector p
    return (toList v)

{- | Experimental implementation of 'luPacked'
     for any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.

>>> let m = ident 5 + (5><5) [0..] :: Matrix (Z ./. 17)
(5><5)
 [  1,  1,  2,  3,  4
 ,  5,  7,  7,  8,  9
 , 10, 11, 13, 13, 14
 , 15, 16,  0,  2,  2
 ,  3,  4,  5,  6,  8 ]

>>> let (l,u,p,s) = luFact $ luPacked' m
>>> l
(5><5)
 [  1,  0, 0,  0, 0
 ,  6,  1, 0,  0, 0
 , 12,  7, 1,  0, 0
 ,  7, 10, 7,  1, 0
 ,  8,  2, 6, 11, 1 ]
>>> u
(5><5)
 [ 15, 16,  0,  2,  2
 ,  0, 13,  7, 13, 14
 ,  0,  0, 15,  0, 11
 ,  0,  0,  0, 15, 15
 ,  0,  0,  0,  0,  1 ]

-}
luPacked' x = LU m p
  where
    (m,p) = mutable (luST (magnit 0)) x

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

scalS a (Slice x r0 c0 nr nc) = rowOper (SCAL a (RowRange r0 (r0+nr-1)) (ColRange c0 (c0+nc-1))) x

view x k r = do
    d <- readMatrix x k k
    let rr = r-1-k
        o  = if k < r-1 then 1 else 0
        s = Slice x (k+1) (k+1) rr rr
        u = Slice x k     (k+1) o  rr
        l = Slice x (k+1) k     rr o
    return (d,u,l,s)

withVec r f = \s x -> do
    p <- newUndefinedVector r
    _ <- f s x p
    v <- unsafeFreezeVector p
    return v

luPacked'' m = (id *** toList) (mutable (withVec (rows m) lu2) m)
  where
    lu2 (r,_) x p = do
        forM_ [0..r-1] $ \k -> do
            pivot x p k
            (d,u,l,s) <- view x k r
            when (magnit 0 d) $ do
                scalS (recip d) l
                gemmm 1 s (-1) l u

    pivot x p k = do
        j <- maxIndex . abs . flatten <$> extractMatrix x (FromRow k) (Col k)
        writeVector p k (j+k)
        swap k (k+j)
      where
        swap i j = rowOper (SWAP i j AllCols) x

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

rowRange m = [0..rows m -1]

at k = Pos (idxs[k])

backSust' lup rhs = foldl' f (rhs?[]) (reverse ls)
  where
    ls  = [ (d k , u k , b k) | k <- rowRange lup ]
      where
        d k = lup ?? (at k, at k)
        u k = lup ?? (at k, Drop (k+1))
        b k = rhs ?? (at k, All)

    f x (d,u,b) = (b - u<>x) / d
                       ===
                        x

forwSust' lup rhs = foldl' f (rhs?[]) ls
  where
    ls  = [ (l k , b k) | k <- rowRange lup ]
      where
        l k = lup ?? (at k, Take k)
        b k = rhs ?? (at k, All)

    f x (l,b) =     x
                   ===
                (b - l<>x)

luSolve'' (LU lup p) b = backSust' lup (forwSust' lup pb)
  where
    pb = b ?? (Pos (fixPerm' p), All)

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

forwSust lup rhs = fst $ mutable f rhs
  where
    f (r,c) x = do
        l <- unsafeThawMatrix lup
        let go k = gemmm 1 (Slice x k 0 1 c) (-1) (Slice l k 0 1 k) (Slice x 0 0 k c)
        mapM_ go [0..r-1]

backSust lup rhs = fst $ mutable f rhs
  where
    f (r,c) m = do
        l <- unsafeThawMatrix lup
        let d k = recip (lup `atIndex` (k,k))
            u k = Slice l k (k+1) 1 (r-1-k)
            b k = Slice m k 0 1 c
            x k = Slice m (k+1) 0 (r-1-k) c
            scal k = rowOper (SCAL (d k) (Row k) AllCols) m

            go k = gemmm 1 (b k) (-1) (u k) (x k) >> scal k
        mapM_ go [r-1,r-2..0]

{- | Experimental implementation of 'luSolve' for any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.

>>> let a = (2><2) [1,2,3,5] :: Matrix (Z ./. 13)
(2><2)
 [ 1, 2
 , 3, 5 ]
>>> b
(2><3)
 [ 5, 1, 3
 , 8, 6, 3 ]

>>> luSolve' (luPacked' a) b
(2><3)
 [ 4,  7, 4
 , 7, 10, 6 ]

-}
luSolve' (LU lup p) b = backSust lup (forwSust lup pb)
  where
    pb = b ?? (Pos (fixPerm' p), All)

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

data MatrixView t b
    = Elem t
    | Block b b b b
  deriving Show

viewBlock' r c m
    | (rt,ct) == (1,1) = Elem (atM' m 0 0)
    | otherwise        = Block m11 m12 m21 m22
  where
    (rt,ct) = size m
    m11 = subm (0,0) (r,c)       m
    m12 = subm (0,c) (r,ct-c)    m
    m21 = subm (r,0) (rt-r,c)    m
    m22 = subm (r,c) (rt-r,ct-c) m
    subm = subMatrix

viewBlock m = viewBlock' n n m
  where
    n = rows m `div` 2

invershur (viewBlock -> Block a b c d) = fromBlocks [[a',b'],[c',d']]
  where
    r1 = invershur a
    r2 = c <> r1
    r3 = r1 <> b
    r4 = c <> r3
    r5 = r4-d
    r6 = invershur r5
    b' = r3 <> r6
    c' = r6 <> r2
    r7 = r3 <> c'
    a' = r1-r7
    d' = -r6

invershur x = recip x

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

instance Testable (Matrix I) where
   checkT _ = test

test :: (Bool, IO())
test = (and ok, return ())
  where
    m  = (3><4) [1..12] :: Matrix I
    r  = (2><3) [1,2,3,4,3,2]
    c  = (3><2) [0,4,4,1,2,3]
    p  = (9><10) [0..89] :: Matrix I
    ep = (2><3) [10,24,32,44,31,23]
    md = fromInt m      :: Matrix Float
    ok = [ tr m <> m == toInt (tr md <> md)
         , m <> tr m == toInt (md <> tr md)
         , m ?? (Take 2, Take 3) == remap (asColumn (range 2)) (asRow (range 3)) m
         , remap r (tr c) p == ep
         , tr p ?? (PosCyc (idxs[-5,13]), Pos (idxs[3,7,1])) == (2><3) [35,75,15,33,73,13]
         ]

max_stars_repo_name: schnecki/hmatrix-float

max_issues_repo_name: schnecki/hmatrix-float

max_forks_repo_name: schnecki/hmatrix-float


text:
{-# LANGUAGE LambdaCase #-}

module Statistics.Quantile.Approximate.Sampling (
    sampling
  , samplingJackknife
) where

import Control.Applicative
import Control.Monad.IO.Class

import qualified Data.List.Ordered as Bag
import Data.Monoid
import Data.Vector.Unboxed (Vector)
import qualified Data.Vector.Unboxed as V

import Pipes
import qualified Pipes.Prelude as P

import qualified Statistics.Quantile as Q
import Statistics.Resampling
import Statistics.Types
import Statistics.Sample

import Statistics.Quantile.Types
import Statistics.Quantile.Util

import System.Random.MWC

data Jackknife = Jackknifed | NotJackknifed
  deriving (Eq, Show)

sampling :: Int
         -> Int
         -> Selector IO
sampling k n = Selector $ sampling' NotJackknifed k n

samplingJackknife :: Int
                  -> Int
                  -> Selector IO
samplingJackknife k n = Selector $ sampling' Jackknifed k n

-- FIXME: this breaks at 'select' with low probability.
-- FIXME: fold with something stricter than a tuple
sampling' :: Jackknife
          -> Int
          -> Int
          -> Quantile
          -> Stream IO
          -> IO Double
sampling' j k n q (Stream src) = withSystemRandom $ \gen -> do
  P.foldM (sample gen) (pure ([],0)) (finalize j) src
  where choose :: GenIO -> IO Bool
        choose g = choose' <$> uniformR (0.0, invP) g

        choose' :: Double -> Bool
        choose' x
          | x < 1.0   = True
          | otherwise = False

        -- 1/P(choose)
        invP :: Double
        invP = (fromIntegral n) / (fromIntegral k)

        sample :: GenIO -> ([Double], Int) -> Double -> IO ([Double], Int)
        sample g (vs,c) v = choose g >>= \case
          True -> let vs' = Bag.insertBag v vs in
                  pure (vs', c+1)
          False -> pure (vs,c)

        finalize NotJackknifed (vs, c) = pure $ vs !! (qidx q c)
        finalize Jackknifed (vs, _) = pure . estimateQ median $ jackknife (Function $ estimateQ q) $ V.fromList vs

estimateQ :: Quantile -> (Vector Double) -> Double
estimateQ (Quantile k n) = Q.weightedAvg k n

max_stars_repo_name: fractalcat/slides

max_issues_repo_name: fractalcat/slides

max_forks_repo_name: fractalcat/slides


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

import           Numeric.LinearAlgebra

import           Grenade
import           Grenade.Layers.Internal.Convolution
import           Grenade.Layers.Internal.Pooling

import           Criterion.Main

main :: IO ()
main = do
  putStrLn $ "Benchmarking with type: " ++ nameF
  x :: S ('D2 60 60) <- randomOfShape
  y :: S ('D3 60 60 1) <- randomOfShape
  defaultMain
    [ bgroup
        "im2col"
        [ bench "im2col 28x28"     $ whnf (im2col 5 5 1 1) ((28 >< 28) [1 ..])
        , bench "im2col 100x100"   $ whnf (im2col 10 10 1 1) ((100 >< 100) [1 ..])
        , bench "im2col 3x416x416" $ whnf (vid2col 3 3 1 1 416 416) (((3 * 416) >< 416) [1 ..])
        ]
    , bgroup
        "col2im"
        [ bench "col2im 28x28"     $ whnf (col2im 5 5 1 1 28 28) ((576 >< 25) [1 ..])
        , bench "col2im 100x100"   $ whnf (col2im 10 10 1 1 100 100) ((8281 >< 100) [1 ..])
        , bench "col2im 3x416x416" $ whnf (col2vid 1 1 1 1 414 414) (((414 * 414) >< 27) [1 ..])
        ]
    , bgroup
        "poolfw"
        [ bench "poolforwards 3x4" $ whnf (poolForward 1 3 4 2 2 1 1) ((3 >< 4) [1 ..])
        , bench "poolforwards 28x28" $ whnf (poolForward 1 28 28 5 5 1 1) ((28 >< 28) [1 ..])
        , bench "poolforwards 100x100" $ whnf (poolForward 1 100 100 10 10 1 1) ((100 >< 100) [1 ..])
        ]
    , bgroup
        "poolbw"
        [ bench "poolbackwards 3x4" $ whnf (poolBackward 1 3 4 2 2 1 1 ((3 >< 4) [1 ..])) ((2 >< 3) [1 ..])
        , bench "poolbackwards 28x28" $ whnf (poolBackward 1 28 28 5 5 1 1 ((28 >< 28) [1 ..])) ((24 >< 24) [1 ..])
        , bench "poolbackwards 100x100" $ whnf (poolBackward 1 100 100 10 10 1 1 ((100 >< 100) [1 ..])) ((91 >< 91) [1 ..])
        ]
    , bgroup
        "padcrop"
        [ bench "pad 2D 60x60" $ nf (testRun2D Pad) x
        , bench "pad 3D 60x60" $ nf (testRun3D Pad) y
        , bench "crop 2D 60x60" $ nf (testRun2D' Crop) x
        , bench "crop 3D 60x60" $ nf (testRun3D' Crop) y
        ]
    ]
  putStrLn $ "Benchmarked with type: " ++ nameF

testRun2D :: Pad 1 1 1 1 -> S ('D2 60 60) -> S ('D2 62 62)
testRun2D = snd ... runForwards

testRun3D :: Pad 1 1 1 1 -> S ('D3 60 60 1) -> S ('D3 62 62 1)
testRun3D = snd ... runForwards

testRun2D' :: Crop 1 1 1 1 -> S ('D2 60 60) -> S ('D2 58 58)
testRun2D' = snd ... runForwards

testRun3D' :: Crop 1 1 1 1 -> S ('D3 60 60 1) -> S ('D3 58 58 1)
testRun3D' = snd ... runForwards

(...) :: (a -> b) -> (c -> d -> a) -> c -> d -> b
(...) = (.) . (.)

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GADTs                  #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
-- | Minimil linear algebra lib.
module Math.Regression.Simple.LinAlg (
    -- * Operations
    Add (..),
    Eye (..),
    Mult (..),
    Det (..),
    Inv (..),
    -- * Zeros
    zerosLin,
    zerosQuad,
    optimaQuad,
    -- * Two dimensions
    V2 (..),
    M22 (..),
    SM22 (..),
    -- * Three dimensions
    V3 (..),
    M33 (..),
    SM33 (..),
) where

import Control.DeepSeq (NFData (..))
import Data.Complex    (Complex (..))

-------------------------------------------------------------------------------
-- Classes
-------------------------------------------------------------------------------

-- | Addition
class Add a where
    zero :: a
    add  :: a -> a -> a

-- | Identity
class Eye a where
    eye :: a

-- | Multiplication of different things.
class Eye a => Mult a b c | a b -> c where
    mult :: a -> b -> c

-- | Determinant
class Eye a => Det a where
    det :: a -> Double

-- | Inverse
class Det a => Inv a where
    inv :: a -> a

infixl 6 `add`
infixl 7 `mult`

instance Eye Double where
    eye = 1

instance Add Double where
    zero = 0
    add = (+)

instance Mult Double Double Double where
    mult = (*)

instance Det Double where
    det = id

instance Inv Double where
    inv = recip

-------------------------------------------------------------------------------
-- Zeros
-------------------------------------------------------------------------------

-- | Solve linear equation.
--
-- >>> zerosLin (V2 1 2)
-- -2.0
--
zerosLin :: V2 -> Double
zerosLin (V2 a b) = negate (b / a)

-- | Solve quadratic equation.
--
-- >>> zerosQuad (V3 2 0 (-1))
-- Right (-0.7071067811865476,0.7071067811865476)
--
-- >>> zerosQuad (V3 2 0 1)
-- Left ((-0.0) :+ (-0.7071067811865476),(-0.0) :+ 0.7071067811865476)
--
-- Double root is not treated separately:
--
-- >>> zerosQuad (V3 1 0 0)
-- Right (-0.0,0.0)
--
-- >>> zerosQuad (V3 1 (-2) 1)
-- Right (1.0,1.0)
--
zerosQuad :: V3 -> Either (Complex Double, Complex Double) (Double, Double)
zerosQuad (V3 a b c)
    | delta < 0 = Left ((-b/da) :+ (-sqrtNDelta/da), (-b/da) :+ (sqrtNDelta/da))
    | otherwise = Right ((- b - sqrtDelta) / da, (-b + sqrtDelta) / da)
  where
    delta = b*b - 4 * a * c
    sqrtDelta = sqrt delta
    sqrtNDelta = sqrt (- delta)
    da = 2 * a

-- | Find an optima point.
--
-- >>> optimaQuad (V3 1 (-2) 0)
-- 1.0
--
-- compare to
--
-- >>> zerosQuad (V3 1 (-2) 0)
-- Right (0.0,2.0)
--
optimaQuad :: V3 -> Double
optimaQuad (V3 a b _) = zerosLin (V2 (2 * a) b)

-------------------------------------------------------------------------------
-- 2 dimensions
-------------------------------------------------------------------------------

-- | 2d vector. Strict pair of 'Double's.
--
-- Also used to represent linear polynomial: @V2 a b@  \(= a x + b\).
--
data V2 = V2 !Double !Double
  deriving (Eq, Show)

instance NFData V2 where
    rnf V2 {} = ()

instance Add V2 where
    zero = V2 0 0
    add (V2 x y) (V2 x' y') = V2 (x + x') (y + y')
    {-# INLINE zero #-}
    {-# INLINE add #-}

instance Mult Double V2 V2 where
    mult k (V2 x y) = V2 (k * x) (k * y)
    {-# INLINE mult #-}

-- | 2×2 matrix.
data M22 = M22 !Double !Double !Double !Double
  deriving (Eq, Show)

instance NFData M22 where
    rnf M22 {} = ()

instance Add M22 where
    zero = M22 0 0 0 0
    add (M22 a b c d) (M22 a' b' c' d') = M22 (a + a') (b + b') (c + c') (d + d')
    {-# INLINE zero #-}
    {-# INLINE add #-}

instance Eye M22 where
    eye = M22 1 0 0 1
    {-# INLINE eye #-}

instance Det M22 where det = det2
instance Inv M22 where inv = inv2

instance Mult Double M22 M22 where
    mult k (M22 a b c d) = M22 (k * a) (k * b) (k * c) (k * d)
    {-# INLINE mult #-}

instance Mult M22 V2 V2 where
    mult (M22 a b c d) (V2 u v) = V2 (a * u + b * v) (c * u + d * v)
    {-# INLINE mult #-}

-- | >>> M22 1 2 3 4 `mult` eye @M22
-- M22 1.0 2.0 3.0 4.0
instance Mult M22 M22 M22 where
    mult (M22 a b c d) (M22 x y z w) = M22
        (a * x + b * z) (a * y + b * w)
        (c * x + d * z) (c * y + d * w)
    {-# INLINE mult #-}

det2 :: M22 -> Double
det2 (M22 a b c d) = a * d - b * c
{-# INLINE det2 #-}

inv2 :: M22 -> M22
inv2 m@(M22 a b c d) = M22
    (  d / detm) (- b / detm)
    (- c / detm) (  a / detm)
  where
    detm = det2 m
{-# INLINE inv2 #-}

-- | Symmetric 2x2 matrix.
data SM22 = SM22 !Double !Double !Double
  deriving (Eq, Show)

instance NFData SM22 where
    rnf SM22 {} = ()

instance Add SM22 where
    zero = SM22 0 0 0
    add (SM22 a b d) (SM22 a' b' d') = SM22 (a + a') (b + b') (d + d')
    {-# INLINE zero #-}
    {-# INLINE add #-}

instance Eye SM22 where
    eye = SM22 1 0 1
    {-# INLINE eye #-}

instance Det SM22 where det = detS2
instance Inv SM22 where inv = invS2

instance Mult Double SM22 SM22 where
    mult k (SM22 a b d) = SM22 (k * a) (k * b) (k * d)
    {-# INLINE mult #-}

instance Mult SM22 V2 V2 where
    mult (SM22 a b d) (V2 u v) = V2 (a * u + b * v) (b * u + d * v)
    {-# INLINE mult #-}

detS2 :: SM22 -> Double
detS2 (SM22 a b d) = a * d - b * b
{-# INLINE detS2 #-}

invS2 :: SM22 -> SM22
invS2 m@(SM22 a b d) = SM22
    (  d / detm)
    (- b / detm) (  a / detm)
  where
    detm = detS2 m
{-# INLINE invS2 #-}

-------------------------------------------------------------------------------
-- 3 dimensions
-------------------------------------------------------------------------------

-- | 3d vector. Strict triple of 'Double's.
--
-- Also used to represent quadratic polynomial: @V3 a b c@  \(= a x^2 + b x + c\).
data V3 = V3 !Double !Double !Double
  deriving (Eq, Show)

instance NFData V3 where
    rnf V3 {} = ()

instance Add V3 where
    zero = V3 0 0 0
    add (V3 x y z) (V3 x' y' z') = V3 (x + x') (y + y') (z + z')
    {-# INLINE zero #-}
    {-# INLINE add #-}

instance Mult Double V3 V3 where
    mult k (V3 x y z) = V3 (k * x) (k * y) (k * z)
    {-# INLINE mult #-}

-- | 3×3 matrix.
data M33 = M33
    !Double !Double !Double
    !Double !Double !Double
    !Double !Double !Double
  deriving (Eq, Show)

instance NFData M33 where
    rnf M33 {} = ()

instance Add M33 where
    zero = M33 0 0 0 0 0 0 0 0 0

    add (M33 a b c d e f g h i) (M33 a' b' c' d' e' f' g' h' i') = M33
        (a + a') (b + b') (c + c')
        (d + d') (e + e') (f + f')
        (g + g') (h + h') (i + i')
    {-# INLINE zero #-}
    {-# INLINE add #-}

instance Eye M33 where
    eye = M33 1 0 0
              0 1 0
              0 0 1
    {-# INLINE eye #-}

instance Det M33 where det = det3
instance Inv M33 where inv = inv3

instance Mult Double M33 M33 where
    mult k (M33 a b c d e f g h i) = M33
        (k * a) (k * b) (k * c)
        (k * d) (k * e) (k * f)
        (k * g) (k * h) (k * i)
    {-# INLINE mult #-}

instance Mult M33 V3 V3 where
    mult (M33 a b c
           d e f
           g h i) (V3 u v w) = V3
        (a * u + b * v + c * w)
        (d * u + e * v + f * w)
        (g * u + h * v + i * w)
    {-# INLINE mult #-}

-- TODO: instance Mult M33 M33 M33 where

det3 :: M33 -> Double
det3 (M33 a b c
          d e f
          g h i)
    = a * (e*i-f*h) - d * (b*i-c*h) + g * (b*f-c*e)
{-# INLINE det3 #-}

inv3 :: M33 -> M33
inv3 m@(M33 a b c
            d e f
            g h i)
    = M33 a' b' c'
          d' e' f'
          g' h' i'
  where
    a' = cofactor e f h i / detm
    b' = cofactor c b i h / detm
    c' = cofactor b c e f / detm
    d' = cofactor f d i g / detm
    e' = cofactor a c g i / detm
    f' = cofactor c a f d / detm
    g' = cofactor d e g h / detm
    h' = cofactor b a h g / detm
    i' = cofactor a b d e / detm
    cofactor q r s t = det2 (M22 q r s t)
    detm = det3 m
{-# INLINE inv3 #-}

-- | Symmetric 3×3 matrix.
data SM33 = SM33
    !Double
    !Double !Double
    !Double !Double !Double
  deriving (Eq, Show)

instance NFData SM33 where
    rnf SM33 {} = ()

instance Add SM33 where
    zero = SM33 0 0 0 0 0 0

    add (SM33 a d e g h i) (SM33 a' d' e' g' h' i') = SM33
        (a + a')
        (d + d') (e + e')
        (g + g') (h + h') (i + i')
    {-# INLINE zero #-}
    {-# INLINE add #-}

instance Eye SM33 where
    eye = SM33 1
               0 1
               0 0 1
    {-# INLINE eye #-}

instance Det SM33 where det = detS3
instance Inv SM33 where inv = invS3

instance Mult Double SM33 SM33 where
    mult k (SM33 a d e g h i) = SM33
        (k * a)
        (k * d) (k * e)
        (k * g) (k * h) (k * i)
    {-# INLINE mult #-}

instance Mult SM33 V3 V3 where
    mult (SM33 a
               d e
               g h i) (V3 u v w) = V3
        (a * u + d * v + g * w)
        (d * u + e * v + h * w)
        (g * u + h * v + i * w)
    {-# INLINE mult #-}

detS3 :: SM33 -> Double
detS3 (SM33 a
            d e
            g h i)
    = a * (e*i-h*h) - d * (d*i-g*h) + g * (d*h-g*e)
{-# INLINE detS3 #-}

invS3 :: SM33 -> SM33
invS3 m@(SM33 a
              d e
              g h i)
     = SM33 a'
            d' e'
            g' h' i'
  where
    a' = cofactor e h h i / detm
    d' = cofactor h d i g / detm
    e' = cofactor a g g i / detm
    g' = cofactor d e g h / detm
    h' = cofactor d a h g / detm
    i' = cofactor a d d e / detm
    cofactor q r s t = det2 (M22 q r s t)
    detm = detS3 m
{-# INLINE invS3 #-}

max_stars_repo_name: phadej/regression-simple

max_issues_repo_name: phadej/regression-simple

max_forks_repo_name: phadej/regression-simple


text:
{-# LANGUAGE CPP, ForeignFunctionInterface #-}
-----------------------------------------------------------------------------
{- |
Module      :  Numeric.GSL.Polynomials
Copyright   :  (c) Alberto Ruiz 2006
License     :  GPL-style

Maintainer  :  Alberto Ruiz (aruiz at um dot es)
Stability   :  provisional
Portability :  uses ffi

Polynomials.

<http://www.gnu.org/software/gsl/manual/html_node/General-Polynomial-Equations.html#General-Polynomial-Equations>

-}
-----------------------------------------------------------------------------
module Numeric.GSL.Polynomials (
    polySolve
) where

import Data.Packed.Internal
import Data.Complex
import System.IO.Unsafe (unsafePerformIO)

#if __GLASGOW_HASKELL__ >= 704
import Foreign.C.Types (CInt(..))
#endif

{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/. For example,
     the three solutions of x^3 + 8 = 0

@\> polySolve [8,0,0,1]
[(-1.9999999999999998) :+ 0.0,
 1.0 :+ 1.732050807568877,
 1.0 :+ (-1.732050807568877)]@

The example in the GSL manual: To find the roots of x^5 -1 = 0:

@\> polySolve [-1, 0, 0, 0, 0, 1]
[(-0.8090169943749475) :+ 0.5877852522924731,
(-0.8090169943749475) :+ (-0.5877852522924731),
0.30901699437494734 :+ 0.9510565162951536,
0.30901699437494734 :+ (-0.9510565162951536),
1.0 :+ 0.0]@

-}
polySolve :: [Double] -> [Complex Double]
polySolve = toList . polySolve' . fromList

polySolve' :: Vector Double -> Vector (Complex Double)
polySolve' v | dim v > 1 = unsafePerformIO $ do
    r <- createVector (dim v-1)
    app2 c_polySolve vec v vec r "polySolve"
    return r
             | otherwise = error "polySolve on a polynomial of degree zero"

foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TVCV

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
module Mandelbrot where

import Graphics.UI.GLUT
import Control.Monad
import Data.Int
import Data.Complex
import System.Exit ( exitWith, ExitCode(ExitSuccess) )
import Data.IORef
import qualified Data.ByteString as B
import           Graphics.Rendering.OpenGL.Capture (capturePPM)
import           Graphics.Image
import           System.Directory (removeFile)
import Data.List.Split (splitOn)

convert :: FilePath -> FilePath -> Bool -> IO ()
convert input output remove = do
  ppm <- readImageRGBA VU input
  writeImage output ppm
  when remove $ removeFile input

iterations = 40

x // y = fromIntegral x / fromIntegral y

-- Divides [a] into [[a], [a], ...] with each sublist of length n,
-- except the last sublist which has length <= n.
chunkify n [] = []
chunkify n xs = let (xs', rest) = splitAt n xs
                in xs' : chunkify n rest

-- Converts a coordinate in screen space to a vertex.
pix2vert (Size w h) (x, y) = Vertex2 ((3 // w * fromIntegral x) - 2.0)
                             ((2 // h * fromIntegral y) - 1.0)

-- List of all of the vertices that represent screen pixels.
vertices :: IO [Vertex2 GLfloat]
vertices = get windowSize >>= \(Size w h) ->
           return $ [pix2vert (Size w h) (x, y) | x <- [0..w-1], y <- [0..h-1]]

-- Gets the color for a number of iterations.
color3 r g b = Color3 r g b
getcolor :: Int -> Int -> Color3 Float
getcolor iterations iter
  | iter == iterations = color3 0.1 0.1 0.1
  | otherwise          = color3 (amt*0.5) (amt * 0.5) (amt*0.5)
  where amt = iter // iterations

-- Returns the number of iterations <= the maximum iterations of the
-- Mandelbrot set at the given vertex.
mandel numIterations horizontalShift verticalShift zoom (Vertex2 r i) =
  length . takeWhile (\z -> Data.Complex.magnitude z <= 2) .
  take numIterations $
  iterate (\z -> z^2 + (((r / zoom) + horizontalShift) :+ ((i / zoom) + verticalShift))) 0

-- plots one point.
drawVert iter horShift verShift zoom v = do
  horShift' <- get horShift
  verShift' <- get verShift
  zoom' <- get zoom
  color . (getcolor iter) $ mandel iter horShift' verShift' zoom' v
  vertex v

-- draws the whole fractal
display :: IORef Int -> IORef GLfloat -> IORef GLfloat -> IORef GLfloat ->  IO ()
display numIter horShift verShift zoom = do
  -- clear [ ColorBuffer ]
  displayCallback $= (vertices >>= (display' numIter horShift verShift zoom) . (\x -> [x]))
  get currentWindow >>= postRedisplay

-- draws all the vertices in slices (to update the display while drawing).
display' :: IORef Int -> IORef GLfloat -> IORef GLfloat -> IORef GLfloat -> [[Vertex2 GLfloat]] -> IO ()
display' numIter horShift verShift zoom chunks = do
  iter <- get numIter
  mapM_ (render iter) chunks
  displayCallback $= display numIter horShift verShift zoom
  -- exitWith ExitSuccess
  where render iter' points = do
          renderPrimitive Points (mapM_ (drawVert iter' horShift verShift zoom) points)
          flush

-- keyboard :: Key -> KeyState -> p1 -> p2 -> IO ()
-- keyboard :: IORef GLfloat -> IORef GLfloat -> Key -> KeyState -> p3 -> p4 -> IO ()
keyboard numIter horShift verShift zoom (Char keyPressed) Down _ _ = case keyPressed of
   ' '   -> postRedisplay Nothing
   '\27' -> exitWith ExitSuccess
   '\100' -> do
     val <- shiftAmount zoom
     horShift $~! (+ val)
     displayCallback $= display numIter horShift verShift zoom
   '\97' -> do
     val <- shiftAmount zoom
     horShift $~! (+ (-val))
     displayCallback $= display numIter horShift verShift zoom
   '\115' -> do
     val <- shiftAmount zoom
     verShift $~! (+ (-val))
     displayCallback $= display numIter horShift verShift zoom
   '\119' -> do
     val <- shiftAmount zoom
     verShift $~! (+ val)
     displayCallback $= display numIter horShift verShift zoom
   '\61' -> do
     val <- get zoom
     zoom $~! (+ val)
     displayCallback $= display numIter horShift verShift zoom
   '\45' -> do
     val <- get zoom
     zoom $~! (+ (-(val / 2)))
     displayCallback $= display numIter horShift verShift zoom
   '\49' -> do
     resetIter numIter
     displayCallback $= display numIter horShift verShift zoom
   '\50' -> do
     resetIter numIter
     numIter $~! (+ (20))
     displayCallback $= display numIter horShift verShift zoom
   '\51' -> do
     resetIter numIter
     numIter $~! (+ (100))
     displayCallback $= display numIter horShift verShift zoom
   '\52' -> do
     resetIter numIter
     numIter $~! (+ (400))
     displayCallback $= display numIter horShift verShift zoom
   '\53' -> do
     resetIter numIter
     numIter $~! (+ (1000))
     displayCallback $= display numIter horShift verShift zoom

   'c' -> do
     horShift' <- get horShift
     verShift' <- get verShift
     zoom' <- get zoom
     let name = show horShift' ++ "," ++
                show verShift' ++ "," ++
                show zoom'

     putStrLn ""
     putStrLn ""
     putStrLn "Writing out file:"
     putStrLn $ "images/" ++ name ++ ".png"
     putStrLn ""
     let ppm = "images/" ++ name ++ ".ppm"
         png = "images/" ++ name ++ ".png"
     (>>=) capturePPM (B.writeFile ppm)
     Mandelbrot.convert ppm png True

   _     -> return ()
keyboard _ _ _ _ _ _ _ _ = return ()

shiftAmount zoom = do
  val <- get zoom
  pure (0.1 / val)

resetIter numIter = do
  val <- get numIter
  numIter $~! (+ (iters + (-val)))

parseFile' (file) = return
  where [x, y, z] = splitOn "," $ (splitOn ".png" file) !! 0
        return = ((read x :: GLfloat), (read y :: GLfloat), (read z :: GLfloat))

parseFile :: [String] -> (GLfloat, GLfloat, GLfloat)
parseFile [justOne] = parseFile' $ (splitOn "images/" justOne) !! 1
parseFile _ = (0.0, 0.0, 1.0)

iters = 20
mandelMain = do
  (_, file) <- getArgsAndInitialize
  mapM putStrLn file
  let (initHor, initVer, initZoom) = parseFile file

  initialDisplayMode $= [ SingleBuffered, RGBMode ]
  initialWindowSize $= Size 1200 1024
  initialWindowPosition $= Position 100 100
  createWindow "Mandelbrot"
  clearColor $= Color4 0 0 0 0
  matrixMode $= Projection
  loadIdentity
  ortho (-2) 1 (-1) 1 (-1) 1
  -- fullScreen
  numIter <- newIORef iters
  horShift <- newIORef initHor
  verShift <- newIORef initVer
  zoom <- newIORef initZoom
  keyboardMouseCallback $= Just (keyboard numIter horShift verShift zoom)
  displayCallback $= (display numIter horShift verShift zoom)
  idleCallback $= Just (idle numIter)
  mainLoop

idle :: IORef Int -> IdleCallback
idle numIter = do
  val <- get numIter
  let newVal = if val < 50 then 5 else 100
  numIter $~! (+ 0)
  postRedisplay Nothing

max_stars_repo_name: jordangedney/mandelbrot-explorer

max_issues_repo_name: jordangedney/mandelbrot-explorer

max_forks_repo_name: jordangedney/mandelbrot-explorer