text:
{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE NoMonomorphismRestriction #-}

-- :set -fobject-code in GHCi
module Picture where

import Foreign hiding (rotate)
import Foreign.C.Types

import Control.Applicative

import Data.Complex

default (Int, Float, Integer)

-- Strict
data Color =
  Color
    { r, g, b, a :: !Float
    }
  deriving (Show)

data CColor =
  CColor
    { rC, gC, bC :: !CUChar
    }
  deriving (Show)

instance Storable CColor where
  alignment _ = 1
  sizeOf _ = 3
  peek ptr =
    CColor <$> peekByteOff ptr 0 <*> peekByteOff ptr 1 <*> peekByteOff ptr 2
  poke ptr (CColor d c i) = do
    pokeByteOff ptr 0 d
    pokeByteOff ptr 1 c
    pokeByteOff ptr 2 i

type Point = (Float, Float)

type Image a = Point -> a

type Region = Image Bool

type Filter c = Image c -> Image c

type FilterC = Filter Color

type ImageC = Image Color

renderRegion :: Region -> Image Color
renderRegion r p =
  if r p
    then black
    else white

cOver :: Color -> Color -> Color
cOver (Color r1 g1 b1 a1) (Color r2 g2 b2 a2) =
  Color (h r1 r2) (h g1 g2) (h b1 b2) (h a1 a2)
  where
    h x1 x2 = x1 + (1 - a1) * x2

bilerpBRBW = bilerpC black red blue white

bilerpC ll lr ul ur (wx, wy) = lerpC wy (lerpC wx ll lr) (lerpC wx ul ur)

over :: ImageC -> ImageC -> ImageC
top `over` bot = \p -> top p `cOver` bot p

instance Semigroup Color where
  (<>) = cOver

instance Monoid Color where
  mempty = invisible

invisible = Color 0 0 0 0

black = Color 0 0 0 1

white = Color 1 1 1 1

red = Color 1 0 0 1

green = Color 0 1 0 1

blue = Color 0 0 1 1

yellow = Color 1 1 0 1

cond :: Image Bool -> Image c -> Image c -> Image c
cond =
  liftA3
    (\a b c ->
       if a
         then b
         else c)

lerpI :: Image Frac -> ImageC -> ImageC -> ImageC
lerpI = liftA3 lerpC

empty = const invisible

whiteI = const white

blackI = const black

redI = const red

blueI = const blue

greenI = const green

yellowI = const yellow

blackWhiteIm reg = cond reg blackI whiteI

blueYellowIm reg = cond reg blueI yellowI

crop :: Region -> FilterC
crop reg im = cond reg im mempty

type Warp = Point -> Point

type Vector = (Float, Float)

translateP :: Vector -> Warp
translateP (dx, dy) (x, y) = (x + dx, y + dy)

scaleP :: Vector -> Warp
scaleP (sx, sy) (x, y) = (sx * x, sy * y)

uscaleP :: Float -> Warp
uscaleP s = scaleP (s, s)

rotateP :: Float -> Warp
rotateP θ (x, y) = (x * cos θ - y * sin θ, y * cos θ + x * sin θ)

udisk :: Region
udisk p = distO p < 1

invWarp warp im = im . warp

translate :: Vector -> Filter c
translate (dx, dy) = invWarp (translateP (-dx, -dy))

scale :: Vector -> Filter c
scale (sx, sy) = invWarp (scaleP (1 / sx, 1 / sy))

-- Laws
-- uscale x . uscale y = uscale (x * y)
uscale :: Float -> Filter c
uscale s = invWarp (uscaleP (1 / s))

rotate :: Float -> Filter c
rotate θ = invWarp (rotateP (-θ))

swirlP :: Float -> Warp
swirlP r p = rotateP (distO p * (2 * pi / r)) p

swirl :: Float -> Filter c
swirl r = invWarp (swirlP (-r))

square n (x, y) = abs x <= n / 2 && abs y <= n / 2

-- Figure to visually check orientation
disks =
  renderRegion $
  f <$> udisk <*> translate (2, 2) udisk <*> translate (3, 2) udisk <*> vstrip
  where
    f a b c d = a || b || c || d

-- Figures from "Functional Images" by Conal Elliot.
fig1 = renderRegion $ vstrip

fig2 = renderRegion $ checker

fig3 = renderRegion $ altRings

fig4 = renderRegion $ polarChecker 10

fig5 = renderRegion $ uscale 0.05 gasket

fig6 = renderFrac $ wavDist

fig7 = bilerpBRBW

fig8 = lerpI wavDist (blackWhiteIm (polarChecker 10)) (blueYellowIm checker)

fig9 = ybRings

fig10 = renderRegion udisk

fig11 = renderRegion $ swirl 1 vstrip

fig13 t = renderRegion (swirlingXPos t)

fig14 = renderRegion $ annulus 0.5

fig17 = renderRegion $ shiftXor 2.6 altRings

fig18 = renderRegion $ xorgon 8 (7 / 4) altRings

-- Typo in original paper, 0.2 5 should be 0.25
fig19 = crop (wedgeAnnulus 0.25 10) ybRings

fig20 = crop (swirl 2 (wedgeAnnulus 0.25 10)) ybRings

fig22 = translate (-20, -20) tiledBilerp

fig23 = renderRegion $ uscale 3.5 (radInvert checker)

fig24 = rippleRad 8 0.3 ybRings

fig25 = \t -> rippleRad 8 (cos t / 2) ybRings

fig26 = rippleRad 8 0.3 $ cropRad 1 $ ybRings

fig27 = cropRad 1 $ rippleRad 8 0.3 $ ybRings

fig28 = swirl 8 $ rippleRad 5 0.3 $ cropRad 5 $ ybRings

-- Typo in original paper, washer (1/2) (pi/2) 1 tiledBilerp
-- TODO: Fix visual mismatch with paper.
fig32 = washer (1 / 2) (pi / 2) tiledBilerp

-- Convert an image to a static animation.
static :: Image c -> Anim c
static = const

-- The main animation to run.
-- mainAnim = static . renderFrac $ mandel
mainAnim :: Anim Color
mainAnim = fig32

calculate (x, y, t) = adjust (mainAnim t') (x', y')
  where
    duration = 5000 -- Duration of the animation in ms
    range = 2 * pi -- Time range of the animation
    -- Time normalized from 0 to 1
    normalizedTime = fromIntegral (t `mod` duration) / fromIntegral duration
    scaledTime = range * normalizedTime
    -- scaledTime = fromIntegral t / 10000
    (x', y', t') = (fromIntegral x, fromIntegral y, scaledTime)
    adjust = adjustToWindow
    -- Adjust an image to a window, by translating, scaling and flipping
    adjustToWindow :: Filter c
    adjustToWindow =
      translate (screenWidth / 2, screenHeight / 2) .
      uscale 60 . flipY
    flipY p (x, y) = p (x, -y)

(screenWidth, screenHeight) = (640, 480)

foreign export ccall fillPixelBuffer :: Ptr CColor -> CInt -> IO ()
fillPixelBuffer arr t = pokeArray arr (map calc l1)
  where
    l1 = [(x, y) | y <- [0 .. screenHeight - 1], x <- [0 .. screenWidth - 1]]
    calc (x, y) =
      CColor (truncate (255 * r)) (truncate (255 * g)) (truncate (255 * b))
      where
        Color r g b _ = calculate (x, y, t)

checker :: Region
checker (x, y) = even (floor x + floor y)

vstrip :: Region
vstrip (x, y) = abs x <= 1 / 2

overlay :: ImageC -> ImageC -> ImageC
overlay = liftA2 (<>)

-- Sierpinski triangle
gasket :: Region
gasket (x, y) = floor x .|. floor y == (floor x :: Integer)

altRings p = even (floor (distO p))

distO (x, y) = sqrt (x * x + y * y)

type Frac = Float -- in [0, 1]

wavDist :: Image Frac
wavDist p = (1 + cos (pi * distO p)) / 2

ybRings = lerpI wavDist blueI yellowI

renderFrac :: Image Frac -> ImageC
renderFrac f (x, y) = lerpC (f (x, y)) black white

lerpC :: Frac -> Color -> Color -> Color
lerpC w (Color r1 g1 b1 a1) (Color r2 g2 b2 a2) =
  Color (h r1 r2) (h g1 g2) (h b1 b2) (h a1 a2)
  where
    h x1 x2 = w * x1 + (1 - w) * x2

lighten x c = lerpC x c white

darken x c = lerpC x c black

type PolarPoint = (Float, Float)

fromPolar :: Point -> PolarPoint
fromPolar (ρ, θ) = (ρ * cos θ, ρ * sin θ)

toPolar :: PolarPoint -> Point
toPolar (x, y) = (distO (x, y), atan2 y x)

polarChecker :: Int -> Region
polarChecker = polarize checker

polarize :: Image c -> Int -> PolarPoint -> c
polarize pat n = pat . sc . toPolar
  where
    sc (ρ, θ) = (ρ, θ * n' / π)
    π = pi
    n' = fromIntegral n

-- Animations
type Time = Float

type Anim c = Time -> Image c

swirlingXPos :: Anim Bool
swirlingXPos t = swirl (t * t) xPos

xPos :: Region
xPos (x, y) = x > 0

-- Region algebra
intersect :: Region -> Region -> Region
intersect = liftA2 (&&)

union :: Region -> Region -> Region
union = liftA2 (||)

xorR :: Region -> Region -> Region
xorR = liftA2 xor
  where
    xor True b = not b
    xor _ b = b

compR :: Region -> Region
compR = fmap not

universeR :: Region
universeR = const True

emptyR :: Region
emptyR = const False

r \\ r' = r `intersect` compR r'

annulus :: Frac -> Region
annulus inner = udisk \\ uscale inner udisk

radReg :: Int -> Region
radReg n = test . toPolar
  where
    test (r, a) = even (floor (a * fromIntegral n / pi))

wedgeAnnulus :: Frac -> Int -> Region
wedgeAnnulus inner n = annulus inner `intersect` radReg n

shiftXor :: Float -> Region -> Region
shiftXor r reg = reg' r `xorR` reg' (-r)
  where
    reg' d = translate (d, 0) reg

-- Typo in the original paper, missing argument g
xorgon :: Int -> Float -> Region -> Region
xorgon n r g = xorRs (rf <$> [0 .. n - 1])
  where
    rf :: Int -> Region
    rf i = translate (fromPolar (r, a)) g
      where
        a = fromIntegral i * 2 * pi / fromIntegral n

xorRs :: [Region] -> Region
xorRs = foldr xorR emptyR

tileP :: Vector -> Warp
tileP (w, h) (x, y) = (wrap' w x, wrap' h y)

tiledBilerp = about (1 / 2, 1 / 2) (tile (1, 1)) bilerpBRBW

type HyperFilter c = Filter c -> Filter c

about :: Point -> HyperFilter c
about (x, y) filt = translate (x, y) . filt . translate (-x, -y)

wrap :: Float -> Float -> Float
wrap w x = w * fracPart (x / w)
  where
    fracPart t = t - fromIntegral (truncate t)

wrap' :: Float -> Float -> Float
wrap' w x = wrap w (x + w / 2) - w / 2

tile :: Vector -> Filter c
tile size = invWarp (tileP size)

swirlP' r = polarWarp (\(ρ, θ) -> (ρ, θ + ρ * (2 * pi / r)))

polarWarp warp = fromPolar . warp . toPolar

radInvertP :: Warp
radInvertP = polarWarp (\(ρ, θ) -> (1 / ρ, θ))

radInvert :: Filter c
radInvert = invWarp radInvertP

rippleRadP :: Int -> Float -> Warp
rippleRadP n s =
  polarWarp $ (\(ρ, θ) -> (ρ * (1 + s * sin (fromIntegral n * θ)), θ))

rippleRad :: Int -> Float -> Filter c
rippleRad n s = invWarp (rippleRadP n (-s))

cropRad :: Float -> FilterC
cropRad r = crop (uscale r udisk)

wiggleRotateP :: Float -> Float -> Time -> Warp
wiggleRotateP cycles θmax t = polarWarp warp
  where
    warp (r, a) = (r, a + θmax * sin (t + dt))
      where
        dt = 2 * pi * cycles * (r - 1 / 2)

wiggleRotate :: Float -> Float -> Time -> Filter c
wiggleRotate cycles θmax t = invWarp (wiggleRotateP cycles θmax t)

washer :: Float -> Float -> ImageC -> Time -> ImageC
washer cycles θmax im t = cropRad 1 $ wiggleRotate cycles θmax t $ im

-- Constants
maxIter :: Int -- Max iterations
maxIter = 750

width :: Int -- Image width
height :: Int -- Image height
width = 400

height = 400

-- Note: aspect ratio of (minX, minY), (maxX, maxY) must
-- match aspect ratio of (width, height)
minX :: Float -- Min x-coordinate of graph
maxX :: Float -- Max x-coordinate of graph
minY :: Float -- Min y-coordinate of graph
maxY :: Float -- Max y-coordinate of graph
-- For the zoomed in part of the Mandelbrot:
--minX = -0.826341244461360116
--maxX = -0.8026423086165848822
--minY = -0.2167936114403439588
--maxY = -0.193094675595568725
--For a full view of the mandelbrot
minX = -2.5

maxX = 1.5

minY = -2

maxY = 2

-- The actual fractal part. It basically works on a matrix, which we
-- will call M, that represents a grid of points on the
-- graph. Essentially, M[i, j] is (xList[j], yList[i])
xList :: [Float]
yList :: [Float]
xList = [minX,(minX + ((maxX - minX) / (fromIntegral width - 1))) .. maxX]

yList =
  reverse [minY,(minY + ((maxY - minY) / (fromIntegral height - 1))) .. maxY]

etaFraction :: Complex Float -> Float
etaFraction z = (log (log (magnitude z))) / log 2

smoothEta :: Int -> Complex Float -> Float -- Smooth escape time algorithm value
smoothEta iter z = (fromIntegral iter - etaFraction z) / fromIntegral maxIter

-- Gets the color for the point, in range [0, 1]
color :: Int -> Complex Float -> Float
color iter z = smoothEta iter z -- Smooth escape time algorithm (and invert)

-- color iter z = fromIntegral iter / fromIntegral maxIter
interpolate :: Float -> Int -- Adds an interpolation curve for interpolating color
interpolate v = truncate ((v ^ 12) * 255) -- Polynomial curve

-- interpolate v = chr (truncate(v * 255)) -- Linear
-- The actual fractal algorithm!
frac :: Complex Float -> Complex Float -> Int -> Int
frac c@(x :+ y) z iter
  | iter >= maxIter = 255 -- never escaped, return color value of 255
  | let p = sqrt ((x - 0.25) ^ 2 + y ^ 2)
     in x <= p - 2 * p ^ 2 + 0.25 && (x + 1) ^ 2 + y ^ 2 <= 0.625 = 255
  | otherwise =
    let z' = z * z + c
     in if ((realPart z') * (realPart z') + (imagPart z') * (imagPart z')) > 4
          then interpolate (color iter z')
          else frac c z' (iter + 1)

mandel :: Image Frac
mandel (x, y) = fromIntegral (frac (x :+ y) (0 :+ 0) 0) / 255

-- dotGridGradient :: Int -> Int -> Float -> Float -> IO Float
-- dotGridGradient ix iy x y = do
--   (v1, v2) <- do p' <- randomIO :: IO Float
--                  let (_, p) = properFraction p'
--                  pure (p, sqrt (1 - p ^ 2))
--   let c1 = dx * v1
--   let c2 = dy * v2
--   let res = (c1 + c2)
--   -- print res
--   pure res
--   where
--     dx = x - fromIntegral ix
--     dy = y - fromIntegral iy

-- perlin :: Image (IO Frac)
-- perlin (x, y) = do
--   let (x0, sx) = properFraction x
--   let (y0, sy) = properFraction y
--   let x1 = x0 + 1
--   let y1 = y0 + 1
--   n0 <- dotGridGradient x0 y0 x y
--   n1 <- dotGridGradient x1 y0 x y
--   let ix0 = lerp n0 n1 sx
--   n0' <- dotGridGradient x0 y1 x y
--   n1' <- dotGridGradient x1 y1 x y
--   let ix1 = lerp n0' n1' sx
--   let res' = lerp ix0 ix1 sy
--   let res = abs res'
--   pure res
--   where
--     lerp a0 a1 w = (1.0 - w) * a0 + w * a1

max_stars_repo_name: siraben/functional-images

max_issues_repo_name: siraben/functional-images

max_forks_repo_name: siraben/functional-images


text:
{-# LANGUAGE OverloadedLists #-}
{-# LANGUAGE Strict #-}
{-# LANGUAGE BangPatterns#-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE RankNTypes#-}
{-# LANGUAGE DeriveGeneric, DeriveAnyClass #-}
module Tasks.Swarm where

import           Numeric.LinearAlgebra (Vector, scale)
import qualified Data.Vector.Storable as VS
import           Control.Lens
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV
import           Data.Vector.Mutable (STVector)
import           Control.Monad.ST
import           Data.List (nub)
import           Control.Monad ((<=<), foldM, when)
import Debug.Trace
import qualified Data.HashTable.ST.Cuckoo as HT
import qualified Data.HashMap as HM
import Data.Vector.Strategies
import GHC.Generics
import Control.Arrow

import Data.GameObject
import Physics.Mechanics
import Physics.Collision
import Data.PhysicsData

import Utils (norm,normalize)

type CollisionGrid s = HT.HashTable s (Int,Int,Int) [Int]

integrityWeight = 0.5
separationWeight = 2
gravityWeight = 0.1

cellSize = 4

vehSize = 2
minDist = 4

off = map ((.) (round.(/cellSize))) [(+ vehSize),(+ negate vehSize)]

addInPl :: forall s. Vector Float -> STVector s Float -> ST s ()
addInPl v mv = VS.zipWithM_ (\i a -> MV.modify mv (+a) i) [0,1,2] v

data BoidComputation = BoidComputation { _object :: GameObject
                                       , _sumSp :: Vector Float
                                       , _separ :: Vector Float
                                       , _nrNeighb :: {-# UNPACK #-} Int} deriving Show

makeLenses ''BoidComputation

ins :: Int -> Maybe [Int] -> (Maybe [Int], ())
ins ind (Just l) = (Just (ind:l),())
ins ind Nothing  = (Just [ind],())

markObject t ind=
  mapM_ (\v -> HT.mutate t (VS.unsafeIndex v 0,VS.unsafeIndex v 1,VS.unsafeIndex v 2) $ ins ind) .
  VS.mapM (\x -> ($ x) <$> off) .
  view location

createGrid :: V.Vector GameObject ->
               ST s (CollisionGrid s, V.Vector PrimitiveComputation)
createGrid v = do
  ht <- HT.newSized 6000
  vec <- V.imapM (process ht) v
  return (ht, vec)
  where
    process ht ind ob = do
      markObject ht ind ob
      return $ prepareCPrim ob

getCollisionData :: CollisionGrid s ->
                     V.Vector PrimitiveComputation ->
                     ST s (HM.Map (Int, Int) [Contact])
getCollisionData g prims =
  HT.foldM (\a p -> return $ collisionsInCell prims a (second nub p)) HM.empty g

getContacts :: V.Vector GameObject -> HM.Map (Int, Int) [Contact]
getContacts v =
  traceShowId $ runST $ createGrid v >>= uncurry getCollisionData

resolveContacts :: V.Vector GameObject -> V.Vector GameObject
resolveContacts =
  getContacts >>= HM.foldWithKey folder id
  where folder k = (.) . resolveCollision k

collisionsInCell :: V.Vector PrimitiveComputation ->
                      HM.Map (Int, Int) [Contact] ->
                      (a, [Int]) ->
                      HM.Map (Int, Int) [Contact]
collisionsInCell prims m (x, ix1:ix2:l) =
  let newM = insertCollisions prims m ix1 ix2
  in collisionsInCell prims newM (x, ix2:l)
collisionsInCell _ m _ = m

insertCollisions :: V.Vector PrimitiveComputation ->
                      HM.Map (Int, Int) [Contact] ->
                      Int -> Int ->
                      HM.Map (Int, Int) [Contact]
insertCollisions prims m ix1 ix2 =
  let key = (ix1, ix2)
      check = V.unsafeIndex prims
  in if HM.member key m
     then m
     else let val = view contacts $
                    detectCollision (check ix1) (check ix2)
          in if null val
             then m
             else HM.insert key val m

fillGrid :: V.Vector GameObject -> ST s (CollisionGrid s, STVector s BoidComputation)
fillGrid v = do
  ht <- HT.newSized 6000
  vec <- V.thaw =<< V.imapM (process ht) v
  return (ht,vec)
  where
    process g ind ob = do
      markObject g ind ob
      return $ initBoidComputation ob

computeForces :: (GameObject -> GameObject) -> V.Vector GameObject -> V.Vector GameObject
computeForces = computeForcesGrid fillGrid

incCount = over nrNeighb (+1)

initBoidComputation v = BoidComputation v [0,0,0] [0,0,0] 0

computeBoid (BoidComputation o v1 v2 n) =
  if n /= 0 then
    let Just p = o^.physics
        integ  = scale (integrityWeight / fromIntegral n) $ v1 - p^.speed
        sep    = scale separationWeight v2
    in putForceOnMid (integ+sep) o
  else o

computeForcesGrid :: (forall s. V.Vector GameObject -> ST s (CollisionGrid s, STVector s BoidComputation)) -> (GameObject -> GameObject) -> V.Vector GameObject -> V.Vector GameObject
computeForcesGrid s f v =
  V.map (f . computeBoid) $ runST $
             V.freeze <=< uncurry computeForcesMap <=< s $ v

computeForcesMap :: CollisionGrid s -> STVector s BoidComputation -> ST s (STVector s BoidComputation)
computeForcesMap m v = HT.foldM (\x -> computeForcesST x . nub . snd) v m

computeForcesST :: STVector s BoidComputation -> [Int] -> ST s (STVector s BoidComputation)
computeForcesST v (i:l@(_:_)) =
  foldM (computeForcesPair i) v l >>= flip computeForcesST l
computeForcesST v _ = pure v

computeForcesPair
  :: Int
  -> STVector s BoidComputation
  -> Int
  -> ST s (STVector s BoidComputation)
computeForcesPair i1 v i2 = do
  obj1 <- view object <$> MV.read v i1
  obj2 <- view object <$> MV.read v i2
  when (norm (obj1^.location - obj2^.location) < minDist) $ do
    let sepForce = sep obj1 obj2
    MV.modify v (over separ (+ sepForce)) i1
    MV.modify v (over sumSp (+ integ obj2)) i1
    MV.modify v (over separ (subtract sepForce)) i2
    MV.modify v (over sumSp (+ integ obj1)) i2
    MV.modify v incCount i1
    MV.modify v incCount i2
  pure v

sep :: GameObject -> GameObject -> Vector Float
sep obj1 obj2 = obj1^.location - obj2^.location
  -- let v = obj1^.location - obj2^.location
  --     m = norm v
  -- in scale (1/m/(m-1)) v

integ :: GameObject -> Vector Float
integ obj =
  case obj^.physics of
    Just x  -> x^.speed
    Nothing -> [0,0,0]

clamp :: Float -> Vector Float -> Vector Float
clamp m v =
  let n = norm v
  in if n > m then scale (m/n) v else v

addGravityForce :: Vector Float -> GameObject -> GameObject
addGravityForce gravCent obj =
  let direction = gravCent - obj^.location
--       dnorm = norm direction
      force = scale gravityWeight
              -- . normalize $
              direction
  in putForceOnMid force obj
    -- over (physics._Just.acc) (+force) obj

ext: hs

max_stars_repo_name: Antystenes/CPG

max_issues_repo_name: Antystenes/CPG

max_forks_repo_name: Antystenes/CPG


text:
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE InstanceSigs          #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecordWildCards       #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
module Grenade.Layers.FullyConnected (
    FullyConnected (..)
  , FullyConnected' (..)
  , randomFullyConnected
  ) where

import           Control.DeepSeq
import           GHC.Generics                    (Generic)
import           GHC.TypeLits

import           System.Random.MWC               hiding (create)
import           Control.Monad.Primitive         (PrimBase, PrimState)

import           Data.Proxy
import           Data.Serialize

import qualified Numeric.LinearAlgebra           as LA
import           Numeric.LinearAlgebra.Static

import           Grenade.Core

import           Grenade.Layers.Internal.Update

-- | A basic fully connected (or inner product) neural network layer.
data FullyConnected i o = FullyConnected
                        !(FullyConnected' i o)   -- Neuron weights
                        !(FullyConnected' i o)   -- Neuron momentum
                        deriving (Generic, NFData)

data FullyConnected' i o = FullyConnected'
                         !(R o)   -- Bias
                         !(L o i) -- Activations
                        deriving (Generic, NFData)

instance Show (FullyConnected i o) where
  show FullyConnected {} = "FullyConnected"

instance (KnownNat i, KnownNat o) => UpdateLayer (FullyConnected i o) where
  type Gradient (FullyConnected i o) = (FullyConnected' i o)

  runUpdate LearningParameters {..} (FullyConnected (FullyConnected' oldBias oldActivations) (FullyConnected' oldBiasMomentum oldMomentum)) (FullyConnected' biasGradient activationGradient) =
    let (newBias, newBiasMomentum)    = descendVector learningRate learningMomentum learningRegulariser oldBias biasGradient oldBiasMomentum
        (newActivations, newMomentum) = descendMatrix learningRate learningMomentum learningRegulariser oldActivations activationGradient oldMomentum
    in FullyConnected (FullyConnected' newBias newActivations) (FullyConnected' newBiasMomentum newMomentum)

instance (KnownNat i, KnownNat o) => Layer (FullyConnected i o) ('D1 i) ('D1 o) where
  type Tape (FullyConnected i o) ('D1 i) ('D1 o) = R i
  -- Do a matrix vector multiplication and return the result.
  runForwards (FullyConnected (FullyConnected' wB wN) _) (S1D v) = (v, S1D (wB + wN #> v))

  -- Run a backpropogation step for a full connected layer.
  runBackwards (FullyConnected (FullyConnected' _ wN) _) x (S1D dEdy) =
          let wB'  = dEdy
              mm'  = dEdy `outer` x
              -- calcluate derivatives for next step
              dWs  = tr wN #> dEdy
          in  (FullyConnected' wB' mm', S1D dWs)

instance (KnownNat i, KnownNat o) => Serialize (FullyConnected i o) where
  put (FullyConnected (FullyConnected' b w) _) = do
    putListOf put . LA.toList . extract $ b
    putListOf put . LA.toList . LA.flatten . extract $ w

  get = do
      let f  = fromIntegral $ natVal (Proxy :: Proxy i)
      b     <- maybe (fail "Vector of incorrect size") return . create . LA.fromList =<< getListOf get
      k     <- maybe (fail "Vector of incorrect size") return . create . LA.reshape f . LA.fromList =<< getListOf get
      let bm = konst 0
      let mm = konst 0
      return $ FullyConnected (FullyConnected' b k) (FullyConnected' bm mm)

instance (KnownNat i, KnownNat o, KnownNat (i*o)) => RandomLayer (FullyConnected i o) where

  createRandomWith = randomFullyConnected

randomFullyConnected :: forall m i o . (PrimBase m, KnownNat i, KnownNat o, KnownNat (i*o))
                     => WeightInitMethod -> Gen (PrimState m) -> m (FullyConnected i o)
randomFullyConnected m gen = do
  wN <- getRandomMatrix i o m gen
  wB <- getRandomVector i o m gen
  let bm = konst 0
      mm = konst 0
  return $ FullyConnected (FullyConnected' wB wN) (FullyConnected' bm mm)
  where i = natVal (Proxy :: Proxy i)
        o = natVal (Proxy :: Proxy o)

-------------------- GNum instances --------------------

instance (KnownNat i, KnownNat o) => GNum (FullyConnected i o) where
  s |* FullyConnected i o = FullyConnected (s |* i) o
  FullyConnected i o |+ FullyConnected i2 o2 = FullyConnected (0.5 |* (i |+ i2)) (0.5 |* (o |+ o2))
  gFromRational r = FullyConnected (gFromRational r) (gFromRational 0)

instance (KnownNat i, KnownNat o) => GNum (FullyConnected' i o) where
  s |* FullyConnected' i o = FullyConnected' (fromRational s * i) o
  FullyConnected' i o |+ FullyConnected' i2 o2 = FullyConnected' (0.5 * (i + i2)) (0.5 * (o + o2))
  gFromRational r = FullyConnected' (fromRational r) (fromRational r)

max_stars_repo_name: koenigmaximilian/grenade

max_issues_repo_name: koenigmaximilian/grenade

max_forks_repo_name: koenigmaximilian/grenade


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ScopedTypeVariables #-}

module ImageLoader
(
  loadExamples
, loadNonExamples
)

where

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

import Control.Exception          ( Exception
                                  , try )
import Control.Monad              ( sequence )
import System.FilePath            ( FilePath
                                  , (</>) )
import Vision.Image.Type          ( manifestVector )
import Vision.Image.RGB.Type      ( RGB )
import Vision.Image.Grey.Type     ( Grey )
import Vision.Image.Storage.DevIL ( Autodetect (..)
                                  , load )
import Paths_img_neural           ( getDataFileName )

import qualified Data.Vector.Storable       as VS
import qualified Numeric.LinearAlgebra.Data as VD

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

type Samples = [(VD.Vector Double, VD.Vector Double)]

-- Load the vectors which represent an example of the letter A
loadExamples :: Exception e => IO (Either e Samples)
loadExamples = try $ do
    images <- loadImages exampleImageNames
    return $ map (\x -> (x, output)) images
  where output = VD.fromList [1 :: Double]

-- Load the vectors which do not represent an example of the letter A
loadNonExamples :: Exception e => IO (Either e Samples)
loadNonExamples = try $ do
    images <- loadImages nonExampleImageNames
    return $ map (\x -> (x, output)) images
  where output = VD.fromList [0 :: Double]

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

-- Lists all images that are examples of the letter A
exampleImageNames :: [FilePath]
exampleImageNames = map ("A" </>) [
                      "A1.png"
                    , "A2.png"
                    , "A3.png"
                    , "A4.png"
                    , "A5.png"
                    , "A6.png" ]

-- Lists all images that are _not_ examples of the letter A
nonExampleImageNames :: [FilePath]
nonExampleImageNames = map ("not-A" </>) [
                         "B1.png"
                       , "B2.png"
                       , "B3.png"
                       , "B4.png"
                       , "B5.png"
                       , "B6.png"
                       , "B7.png" ]

-- Loads images from the given list
loadImages :: [FilePath] -> IO [VD.Vector Double]
loadImages fp = do
    fnames  <- map' getFileName fp
    vectors <- map' loadFile fnames
    return $ map convertToLinearVector vectors
  where map' f xs   = sequence $ map f xs
        getFileName = getDataFileName . ("res" </>)

-- Loads the given file to an image, converted to a Vector
loadFile :: FilePath -> IO (VS.Vector Double)
loadFile f = do
    dat <- load Autodetect f
    case dat
      of Left  e           -> error  $ show e
         Right (i :: Grey) -> return $ greyManifestToVector i

-- Converts a Grey into a Vector Int
greyManifestToVector :: Grey -> VS.Vector Double
greyManifestToVector m = VS.map (\i -> 1 - (fromIntegral $ round $ (fromIntegral i)/256)) (manifestVector m)

-- Converts a Data/Vector to a LinearAlgebra/Vector
convertToLinearVector :: VS.Storable a => VS.Vector a -> VD.Vector a
convertToLinearVector = VD.fromList . VS.toList


text:
module Echo where

import Data.Complex
import Linear.Metric
import Linear.V3

echo :: (Floating a) => (a -> Complex a) -> V3 a -> V3 a -> a -> Complex a
echo signal carrier obj t = signal (t - 2 * d / c)
  where
    d = distance carrier obj
    c = 299792458

max_stars_repo_name: nihlete/sar-signal

max_issues_repo_name: nihlete/sar-signal

max_forks_repo_name: nihlete/sar-signal


text:
{-# LANGUAGE OverloadedLists #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE Strict #-}
{-# LANGUAGE TupleSections #-}
module Physics.Mechanics where

import           Numeric.LinearAlgebra (Matrix, Vector, scale, (#>), cross, cmap,dot, ident, tr, (<>), norm_2)
import qualified Number.Quaternion as Q
import           Control.Lens hiding ((<.>))
import qualified Data.Vector.Storable as VS
import           Control.Arrow hiding ((<+>))
import           Debug.Trace
import           Data.Maybe (fromMaybe)
import           Data.Function (on)

import Data.PhysicsData
import Data.GameObject
import Data.Mesh
import Physics.Collision
import qualified Data.Vector as V
import Data.List.Lens

import Utils(step,norm)
import Utils.Quaternions

infixr 1 <+>
(<+>) :: Vector Float -> Vector Float -> Vector Float
(<+>) = VS.zipWith (+)

-- (<*>) :: Vector Float -> Vector Float -> Vector Float
-- (<*>) = VS.zipWith (*)

(<.>) :: Vector Float -> Vector Float -> Float
(<.>) = (VS.sum .) . VS.zipWith (*)

addMaybeForce :: Maybe (Vector Float) -> Maybe PhysicsData -> Maybe PhysicsData
addMaybeForce mf mph = over acc . (<+>) <$> mf <*> mph

putForceLoc :: Vector Float -> Vector Float -> GameObject -> GameObject
putForceLoc point force obj =
          let m    = preview (physics._Just.massInv) obj
              lin  = flip scale (scale (force <.> point) force) <$> m
              ang  = cross force point
          in over physics (addMaybeForce lin)
           . over (physics._Just.torque) (<+>ang)
           $ obj

putForce :: Vector Float -> Vector Float -> GameObject -> GameObject
putForce point force obj =
          let p    = point - obj^.location
          in putForceLoc p force obj

putForceOnMid = over (physics._Just.acc) . (+)

putForceOnTip vec obj =
  let quat = obj^.mesh.quaternion
      tip  = rotateWithQuat quat [0,0,-0.1] + obj^.location
  in putForce tip vec obj

putForceOnTipLoc vec obj =
  let quat = obj^.mesh.quaternion
      f    = rotateWithQuat quat vec
  in putForceOnTip f obj

putForceOnSide vec obj =
  let quat = obj^.mesh.quaternion
      tip  = rotateWithQuat quat [1,0,0] <+> obj^.location
  in putForce tip vec obj

putForceOnSideLoc vec obj =
  let quat = obj^.mesh.quaternion
      f    = rotateWithQuat quat vec
  in putForceOnSide f obj

addScaled
  :: Setter GameObject GameObject (VS.Vector Float) (VS.Vector Float)
  -> Lens PhysicsData PhysicsData (VS.Vector Float) (VS.Vector Float)
  -> GameObject -> GameObject
addScaled lens1 lens2 = applyWithPhysics helper
  where helper ph = over lens1 (<+> VS.map (*step) (ph^.lens2))

applySpeed = addScaled location speed

gForce :: Vector Float
gForce = [0,9,0]

gravityF = over (physics._Just.acc) (subtract gForce)

gravCoeff = scale step gForce

getVec :: Lens PhysicsData PhysicsData (VS.Vector Float) (VS.Vector Float) -> GameObject -> VS.Vector Float
getVec l = fromMaybe [0,0,0] . preview (physics._Just.l)

getMat :: Lens PhysicsData PhysicsData (Matrix Float) (Matrix Float) -> GameObject -> Matrix Float
getMat l = fromMaybe (ident 3) . preview (physics._Just.l)

floatOnGround obj =
  let floorNormal = [0,1,0]
      loc       = obj^.location
      d         = dot floorNormal loc + 2
      f         = negate (2-d)^(2 :: Int)
  in if d < 1
  then putForceOnMid (scale f gForce) obj
  else obj

calculateRelPos contact obj = contact^.contactPoint - obj^.location

linearInertia obj =
  fromMaybe 0 $ obj^?physics._Just.massInv

adjustMove relPos p@(l, an) =
  let limit = norm relPos * 0.1
      newAn = signum an * limit
  in if abs an > limit
     then (l + an - newAn,newAn)
     else p

collisionMoves :: Contact -> (GameObject, Vector Float) -> (GameObject, Vector Float) -> ((Float, Vector Float), (Float, Vector Float))
collisionMoves contact p1 p2 =
  let ang1 = angularInertia contact p1
      lin1 = linearInertia $ fst p1
      ang2 = angularInertia contact p2
      lin2 = linearInertia $ fst p2
      coef = contact^.penetration / (ang1 + ang2 + lin1 + lin2)
      (fli1, fan1) = adjustMove (snd p1) (coef*lin1, coef*ang1)
      (fli2, fan2) = adjustMove (snd p2) (coef*lin2, coef*ang2)
      finang1 = getRotationPM contact p1 fan1 ang1
      finang2 = getRotationPM contact p1 fan2 ang2
  in ((fli1, finang1),(fli2, finang2))

getRotationPM :: Contact -> (GameObject, Vector Float) -> Float -> Float -> Vector Float
getRotationPM contact (obj, relPos) ang angIn =
  if angIn == 0
  then [0,0,0]
  else
    scale (ang) .
    scale (1/angIn) $
    getInertiaInv obj #> cross relPos (contact^.contactNormal)

applyMoves :: Contact -> (Float, Vector Float) -> GameObject -> GameObject
applyMoves contact (lin, ang) =
  over (mesh.quaternion) (qvMul ang) .
  over location (+ scale lin (contact^.contactNormal))

calculateDeltaVel :: Contact -> (GameObject, Vector Float) -> Float
calculateDeltaVel contact p@(obj, _) =
  case obj^.physics of
    Just pd ->
      (pd^.massInv +) $
      angularInertia contact p
      -- dot (contact^.contactNormal) .
      -- flip cross relPos .
      -- ((pd^.inertiaInv) #>) .
      -- flip cross (contact^.contactNormal) $ relPos
    Nothing -> 0

angCoeff = 1

angularInertia :: Contact ->  (GameObject, Vector Float) -> Float
angularInertia contact (obj, relPos) =
  case obj^.physics of
    Just pd ->
--      (*angCoeff) .
      dot (contact^.contactNormal) .        -- get component of speed pointing along the normal
      (\angS -> cross angS relPos) .        -- get linear speed of pt
      ((pd^.inertiaInv) #>) $               -- Get ang speed
      cross relPos (contact^.contactNormal) -- for unit impulse
    Nothing -> 0

calculateClosingVel :: Contact -> (GameObject, Vector Float) -> Float
calculateClosingVel contact (obj, relPos) =
  case obj^.physics of
    Just pd ->
      -- (makeContactBasis contact #>) .
      dot (contact^.contactNormal) .
      (pd^.speed +) . cross (pd^.angularS) $ relPos
    Nothing -> 0

calculateImpulse :: Contact -> (GameObject, Vector Float) -> (GameObject, Vector Float) -> Vector Float
calculateImpulse contact o1 o2 =
  let deltaV = (("Delta: "++) . show >>= trace) $
               ((+) `on` calculateDeltaVel contact) o1 o2
      closiV = (("Closing: "++) . show >>= trace) $
               ((+) `on` calculateClosingVel contact) o1 o2
      contV  = negate $ closiV
      restitution = if abs contV < 2
                    then 0.0
                    else 0.0
      impulse = -- if signum contV /= signum (contact^.sign)
                -- then
                  scale (contV * (1.0 + restitution)/deltaV) $ contact^.contactNormal
                -- else [0,0,0]
  in impulse

applyImpulse :: Vector Float -> (GameObject, Vector Float) -> GameObject
applyImpulse impulse (obj, relPos) =
  let vChange = -- impulse
        (("Linear impulse: "++) . show >>= trace) .
        fromMaybe [0,0,0] $
        flip scale impulse <$> (obj^?physics._Just.massInv)
      rotChange = (("Angular impulse: "++) . show >>= trace) .
--        scale angCoeff .
        fromMaybe [0,0,0] $
        (#> cross impulse relPos) <$> (obj^?physics._Just.inertiaInv)
  in over (physics._Just.speed) (+ vChange) .
     over (physics._Just.angularS) (+ rotChange) $ obj

resolveCollision :: (Int, Int) -> [Contact] -> V.Vector GameObject -> V.Vector GameObject
resolveCollision (ix1, ix2) (contact:cs) objs =
  let obj1 = objs V.! ix1
      obj2 = objs V.! ix2
      movement = scale (contact^.penetration*contact^.sign) $
                 contact^.contactNormal
      move = case obj1^.physics of
               Just _ -> over (ix 0._2.location) (+movement)
               Nothing -> over (ix 1._2.location) (subtract movement)
      relPos1 = (("RelPos1: "++) . show >>= trace) $ calculateRelPos contact obj1
      relPos2 = (("RelPos2: "++) . show >>= trace) $ calculateRelPos contact obj2
      impulse = (("Impulse: "++) . show >>= trace) .
                scale (contact^.sign) $
                calculateImpulse contact (obj1, relPos1) (obj2, relPos2)
      (m1, m2) = collisionMoves contact (obj1, relPos1) (obj2, relPos2)
      update i p =
        -- applyMoves contact p .
        applyImpulse i
      updates = move $
                [(ix1, update (impulse) m1 (obj1, relPos1))
                ,(ix2, update (-impulse) m2 (obj2, relPos2))]
      newObjs = V.unsafeUpd objs updates
      check   = V.unsafeIndex newObjs
      -- newColls =
      --   take (length cs) $
      --   detectCollision (prepareCPrim $ check ix1) (prepareCPrim $ check ix2)^.contacts
  in
    resolveCollision (ix1, ix2) cs newObjs
resolveCollision _ [] objs = objs

floorReaction obj =
  let floorHeight = [0,-2,0]
      floorNormal = [0,1,0] :: Vector Float
      vert      = obj^.mesh.vertices
      loc       = obj^.location
      r         = quatToMat3 $ obj^.mesh.quaternion
      m         = (1/) $ fromMaybe 1 $ obj^?physics._Just.massInv
      translate = (+loc) . subtract floorHeight
      collides = filter ((<0) . snd)
               . map (id &&& (dot floorNormal . translate) . (tr r #>))
               $ vert
      minD :: Float
      minD = negate . foldr min 0 . map snd $ collides
      response (point, depth) =
        let
          vrel = dot floorNormal $ getVec speed obj + cross (getVec angularS obj) point
          aobj = dot floorNormal $ cross (i #> cross point floorNormal) point
          i    = getMat inertiaInv obj
          j    = negate (0.8) * vrel/(m+aobj)
          force = scale j floorNormal
          torque = getMat inertiaInv obj #> cross point force
            -- getMat inertiaInv obj #> cross point force
        in (force, torque)
      (f,t) = foldr (\(a,b) (c,d) -> (a+c,b+d)) ([0,0,0],[0,0,0])
            . map response
            $ collides
      n = fromIntegral $ length collides
  in if n /= 0
  then over (physics._Just.angularS) (+ scale (1/n) t)
     . over (physics._Just.speed) (+ scale (1/n) f)
     . over (location) (+ scale minD floorNormal)
     $ obj
  else obj

applyAcc :: GameObject -> GameObject
applyAcc = addScaled (physics._Just.speed) acc

applyForce :: GameObject -> GameObject
applyForce = applyWithPhysics helper
  where helper ph obj =
         let tq    = ph^.torque
             om    = ph^.angularS
             inInv = ph^.inertiaInv
             newOm = om <+> (inInv #> scale step tq)
             quatD = qvMul (scale step newOm) $ obj^.mesh.quaternion
             newQ  = Q.normalize . addQuat quatD $ obj^.mesh.quaternion
             inInD = quatToLInv (ph^.inertiaD) newQ
         in set (physics._Just.torque) [0,0,0]
          . set (physics._Just.acc) [0,0,0]
          . applySpeed
--          . floorReaction
          . applyAcc
--          . floatOnGround
          . set (physics._Just.angularS) newOm
          . set (physics._Just.inertiaInv) inInD
          . set (mesh.quaternion) newQ
          $ obj

applyLinear :: GameObject -> GameObject
applyLinear = set (physics._Just.torque) [0,0,0]
              . set (physics._Just.acc) [0,0,0]
              . applySpeed
              -- . floorReaction
              . applyAcc
              -- . gravityF

applyWithPhysics :: (PhysicsData -> GameObject -> GameObject)
                   -> GameObject -> GameObject
applyWithPhysics f obj = case obj^.physics of
  Just phy -> f phy obj
  Nothing  -> obj

maxspeed = 50

yolo :: Float -> Vector Float -> Vector Float
yolo maxVal vec =
  let vecNorm = norm vec
      coeff   = (vecNorm/maxVal)**(1/2)
  in scale (0 - max 1 coeff) vec

friction = applyWithPhysics helper
  where helper phy =
          let maxspeed  = 1000000
              maxtorque = 100
              coeff1    = scale step . yolo maxspeed $ phy^.speed
              coeff2    = scale step . yolo maxtorque $ phy^.angularS
          in over (physics._Just.angularS) (<+> coeff2)
            . over (physics._Just.speed) (<+> coeff1)

  -- applyWithPhysics helper
  -- where helper phy obj =
  --         let quat   = obj^.mesh.quaternion
  --             rot    = rotateWithQuat quat [0,0,-1]
  --               -- (\(a,b,c) -> [a,b,c]) . imag . flip quatConcat (qInverse quat) $ quatConcat quat (0 +:: (0,0,-1))
  --             force  = 40 -- exp (maxspeed - sNorm) - 1
  --         in over (physics._Just.acc) (<+>cmap (*force) rot) obj


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

module Services (postCompute, postLogin, ApiApp) where

import           Config
import           Control.Monad.IO.Class    (liftIO)
import           Control.Monad.Logger      (runNoLoggingT)
import           Control.Monad.Trans
import           Data.Aeson                hiding (json)
import           Data.Complex
import           Data.IORef
import           Data.Monoid
import           Data.Text                 (pack)
import qualified Data.Text                 as T
import           Data.Time.Clock           (getCurrentTime)
import           Database.Persist.Sqlite   (SqlBackend, createSqlitePool,
                                            runMigration, runSqlPool)
import           Db                        (createToken, findToken, migrateAll,
                                            removeOldTokens, runSQL)
import           Eval                      (eval)
import           GHC.Generics              (Generic)
import           Network.HTTP.Types.Status
import           Text.Read                 (readMaybe)
import           Text.StringRandom         (stringRandomIO)
import           Web.Spock
import           Web.Spock.Config
import Data

type MySession = ()
type MyAppState = Config

type ApiAction ctx a = SpockActionCtx  ctx SqlBackend MySession MyAppState a
type ApiApp ctx = SpockCtxM ctx SqlBackend MySession MyAppState ()

postLogin :: ActionCtxT () (WebStateM SqlBackend MySession MyAppState) b
postLogin = do
  usr <- jsonBody' :: ApiAction () User
  case usr of
    User "admin" "admin" -> do
      currentTime <- liftIO getCurrentTime
      token <- liftIO $ stringRandomIO "[a-zA-Z0-9]{50}"
      runSQL $ createToken token currentTime
      json $ object ["token" .= String token]
    _ -> do setStatus status401
            json $ object ["fault" .= String "invalid login and password"]

postCompute :: ActionCtxT () (WebStateM SqlBackend MySession MyAppState) b
postCompute = do
  (ComputeData token expression) <- jsonBody' :: ApiAction () ComputeData
  currentTime <- liftIO getCurrentTime
  runSQL $ removeOldTokens currentTime -- cleanup old tokens
  tokenLine <- runSQL $ findToken (pack token)
  case tokenLine of
    Just _  -> do
      v <- liftIO $ eval expression
      case v of
        Left _ -> do setStatus status418
                     json $ object ["fault" .= String ( pack $ "can't parse expression [" ++ expression ++ "]")]
        Right a -> case a of
          "True" -> json $ object ["result" .= Bool True]
          "False" -> json $ object ["result" .= Bool False]
          _ -> case readMaybe a of
            Nothing -> do setStatus status418
                          json $ object ["fault" .= String ( pack $ "can't parse expression [" ++ expression ++ "]")]
            Just nv -> json $ object ["result" .= Number nv]
    Nothing -> do setStatus status401
                  json $ object ["fault" .= String "invalid token"]


text:
module SpaceState.City(gotoCity, catapult)
where

import Data.List hiding (concat)
import Data.Maybe
import Data.Foldable
import Control.Monad hiding (mapM_)
import Control.Monad.State as State hiding (mapM_)
import Prelude hiding (catch, concat, mapM_)
import Text.Printf

import Graphics.Rendering.OpenGL as OpenGL
import Graphics.Rendering.FTGL as FTGL
import Graphics.UI.SDL as SDL hiding (flip)

import Statistics
import OpenGLUtils
import Politics
import Entity
import AObject
import Cargo
import Utils
import TextScreen
import Mission
import SpaceState.Game

updateAvailableMission :: AObject -> StateT SpaceState IO ()
updateAvailableMission lc = do
  state <- State.get
  let malleg = colonyOwner lc
  case malleg of
    Nothing -> return ()
    Just alleg -> do
      case possibleMissionType (allegAttitude alleg state) of
        Nothing        -> return ()
        Just m         -> createMission m alleg lc

createMission :: MissionCategory -> String -> AObject -> StateT SpaceState IO ()
createMission Messenger alleg lc = do
  let planetname = aobjName lc
  state <- State.get
  let otherplanets = fmap aobjName $ filter (hasOwner alleg) (toList $ aobjects state)
  when (not (null otherplanets)) $ do
    n <- liftIO $ chooseIO otherplanets
    when (n /= planetname) $
      modify $ modAvailMission $ const $ Just $ MessengerMission n

createMission SecretMessage alleg lc = do
  let planetname = aobjName lc
  state <- State.get
  let otherplanets = map aobjName $ filter (\a ->  hasSomeOwner a && (not . hasOwner alleg) a) (toList $ aobjects state)
  when (not (null otherplanets)) $ do
    n <- liftIO $ chooseIO otherplanets
    when (n /= planetname) $
      modify $ modAvailMission $ const $ Just $ SecretMessageMission n

gotoCity :: AObject -> StateT SpaceState IO ()
gotoCity lc = do
  let planetname = aobjName lc
  state <- State.get
  nmarket <- if planetname == fst (lastmarket state)
               then return $ lastmarket state
               else do
                 m <- liftIO $ randomMarket
                 return (planetname, m)
  modify $ modMarket $ const nmarket
  handleArrival lc
  updateAvailableMission lc
  cityLoop lc
  catapult lc

cityLoop :: AObject -> StateT SpaceState IO ()
cityLoop lc = do
  let planetname = aobjName lc
  state <- State.get
  let f = gamefont state
      alleg = getAllegiance lc
      leave = "Leave " ++ planetname
      cangovernor = isJust $ possibleMissionType (allegAttitude alleg state)
      options = "Market" : "Shipyard" : "Foreign affairs" : (if cangovernor then "Governor" : [leave] else [leave])
  n <- liftIO $ menu 1 (f, Color4 1.0 1.0 1.0 1.0,
                  concat ["Starport on " ++ planetname,
                          if alleg == planetname then "" else "\nThis planet belongs to the country of " ++ alleg ++ "."])
            (map (\s -> (f, Color4 1.0 1.0 0.0 1.0, s)) options)
            (f, Color4 1.0 1.0 0.0 1.0, "=>")
  case n of
    1 -> gotoMarket planetname >> cityLoop lc
    2 -> gotoShipyard >> cityLoop lc
    3 -> liftIO (showForeignAffairs (colonyOwner lc) (gamefont state)) >> cityLoop lc
    4 -> if cangovernor then gotoGovernor lc >> cityLoop lc else return ()
    _ -> return ()

showForeignAffairs :: Maybe String -> Font -> IO ()
showForeignAffairs mst f =
  let strs =
        case mst of
          Nothing -> ["This colony is in anarchy!"]
          Just st -> map relToStr rs
            where rs = relationshipsOf st
                  relToStr (l, v) | v <= (-4) = printf "We despise the state of %s." l
                  relToStr (l, v) | v <= (-1) = printf "We are not very fond of %s." l
                  relToStr (l, v) | v <= 1 = printf "Our relationship with the state of %s is neutral." l
                  relToStr (l, v) | v <= 4 = printf "%s is a friend of ours." l
                  relToStr (l, _) | otherwise = printf "We are allied with the state of %s." l
      drawfunc = makeTextScreen (100, 420)
                                (map (\s -> (f, Color4 1.0 1.0 1.0 1.0, s))
                                     (strs ++ ["\n\nPress any key to continue"]))
                                (return ())
  in pressAnyKeyScreen drawfunc

gotoShipyard :: StateT SpaceState IO ()
gotoShipyard = do
  state <- State.get
  let f = gamefont state
  let damages = startPlHealth - plhealth state
  let cost = damages * 20
  let repairtext =
        if damages == 0
          then "Repair (no damages)"
          else "Repair ship (Cost: " ++ show cost ++ ")"
  n <- liftIO $ menu 1 (f, Color4 1.0 1.0 1.0 1.0, "Shipyard")
            [(f, Color4 1.0 1.0 0.0 1.0, repairtext),
             (f, Color4 1.0 1.0 0.0 1.0, "Exit")]
            (f, Color4 1.0 1.0 0.0 1.0, "=>")
  case n of
    1 -> do
           when (cost > 0 && plcash state >= cost) $ do
             modify $ modPlCash $ (subtract cost)
             modify $ modPlHealth $ const startPlHealth
           gotoShipyard
    _ -> return ()

gotoMarket :: String -> StateT SpaceState IO ()
gotoMarket planetname = do
  state <- State.get
  let market = snd . lastmarket $ state
  (m', cargo', cash', hold') <- liftIO $ execStateT
                                  (tradeScreen ("Market on " ++ planetname)
                                      (gamefont state) (monofont state))
                                  (market, plcargo state, plcash state, plholdspace state)
  modify $ modMarket $ modSnd $ const m'
  modify $ modPlCargo $ const cargo'
  modify $ modPlCash $ const cash'
  modify $ modPlHoldspace $ const hold'

catapult :: AObject -> StateT SpaceState IO ()
catapult lc = do
  state <- State.get
  let objpos = AObject.getPosition lc
      dist = AObject.size lc
      plloc = Entity.position (tri state)
      pldir = OpenGLUtils.normalize (plloc *-* objpos)
  modify $ modTri $ modifyPosition $ const $ objpos *+* (pldir *** (dist + 5))
  modify $ modTri $ modifyVelocity $ const $ pldir *** 0.2
  modify $ modTri $ resetAcceleration
  modify $ modTri $ modifyRotation $ (+180)

gotoGovernor :: AObject -> StateT SpaceState IO ()
gotoGovernor lc = do
  state <- State.get
  let alleg = getAllegiance lc
  case missionFor alleg (plmissions state) of
    Nothing          ->
      case availmission state of
        Nothing  -> noMissionsScreen
        Just mis -> offerMission mis lc
    Just currmission -> missionDisplay (gamefont state) alleg currmission

missionDisplay f alleg (MessengerMission tgt) = do
  missionDisplayGeneric f alleg ["You have an important message that you need to",
                                 "bring to the planet of " ++ tgt ++ "."]
missionDisplay f alleg (SecretMessageMission tgt) = do
  missionDisplayGeneric f alleg ["There's a spy of ours currently residing",
                                 "on the foreign planet of " ++ tgt ++ ".",
                                 "You must bring him a secret message."]

missionDisplayGeneric f alleg msg =
  pressAnyKeyScreen
    (liftIO $ makeTextScreen (100, 500)
        [(f, Color4 1.0 0.2 0.2 1.0,
           intercalate "\n" (["Your current mission for " ++ alleg ++ ":", ""]
                             ++ msg ++ ["", "Press any key to continue"]))]
        (return ()))

handleArrival :: AObject -> StateT SpaceState IO ()
handleArrival lc = do
  state <- State.get
  mapM_ (flip handleArrival' lc) (missions (plmissions state))

noMissionsScreen :: StateT SpaceState IO ()
noMissionsScreen = do
  state <- State.get
  let plname = playername state
  pressAnyKeyScreen
    (liftIO $ makeTextScreen (100, 500)
        [(gamefont state, Color4 1.0 0.2 0.2 1.0,
           intercalate "\n" ["\"Dear " ++ plname ++ ", we currently have",
                             "no tasks for you.\"",
                             "",
                             "",
                             "Press any key to continue"])]
        (return ()))

offerMission :: Mission -> AObject -> StateT SpaceState IO ()
offerMission mis@(MessengerMission tgt) lc = do
  state <- State.get
  let plname = playername state
  let txt = ["\"Dear " ++ plname ++ ", I welcome you to my",
            "residence. I have an important mission that",
            "needs to be taken care of by a trustworthy",
            "adventurer like yourself.",
            "",
            "The mission requires you to deliver an",
            "important message from here to the planet",
            "of " ++ tgt ++ ".",
            "Will you accept?\" (y/n)"]
  offerMissionGeneric mis lc txt

offerMission mis@(SecretMessageMission tgt) lc = do
  state <- State.get
  let plname = playername state
  let txt = ["\"Dear " ++ plname ++ ", I welcome you to my",
            "residence. I have a very important mission that",
            "needs to be taken care of by a competent",
            "adventurer like yourself.",
            "",
            "The mission requires you to deliver a message",
            "to a spy of ours hidden on the foreign planet",
            "of " ++ tgt ++ ".",
            "Will you accept?\" (y/n)"]
  offerMissionGeneric mis lc txt

offerMissionGeneric :: Mission -> AObject -> [String] -> StateT SpaceState IO ()
offerMissionGeneric mis lc msg = do
  state <- State.get
  let alleg = getAllegiance lc
  let txt = intercalate "\n" msg
  c <- pressOneOfScreen
             (liftIO $ makeTextScreen (100, 500)
                         [(gamefont state, Color4 1.0 1.0 1.0 1.0, txt)]
                          (return ()))
             [SDLK_y, SDLK_n, SDLK_ESCAPE]
  modify $ modAvailMission (const Nothing)
  case c of
    SDLK_y -> modify $ modPlMissions $ setMission alleg mis
    _      -> return ()

handleArrival' :: (String, Mission) -> AObject -> StateT SpaceState IO ()
handleArrival' (alleg, MessengerMission tgt) lc = do
  state <- State.get
  let plname = playername state
      planetname = aobjName lc
  checkArrived tgt planetname alleg pointsForMessengerMission
    ["\"Thank you, " ++ plname ++ ", for serving our great",
     "country of " ++ alleg ++ " by delivering this important",
     "message to us.\"",
     "",
     "Mission accomplished!",
     "",
     "",
     "Press Enter to continue"]

handleArrival' (alleg, SecretMessageMission tgt) lc = do
  state <- State.get
  let plname = playername state
      planetname = aobjName lc
  checkArrived tgt planetname alleg pointsForSecretMessage
    ["\"Thank you, " ++ plname ++ ", for serving our great",
     "country of " ++ alleg ++ " by delivering this important",
     "message to us here on " ++ planetname ++ ".\"",
     "",
     "Mission accomplished!",
     "",
     "",
     "Press Enter to continue"]

onActualArrival :: String -> Int -> [String] -> StateT SpaceState IO ()
onActualArrival alleg pts str = do
  state <- State.get
  let txt = intercalate "\n" str
  modify $ modAvailMission (const Nothing)
  modify $ modAllegAttitudes $ modAttitude (+pts) alleg
  modify $ modPlMissions $ removeMission alleg
  pressKeyScreen
    (liftIO $ makeTextScreen (100, 500)
        [(gamefont state, Color4 1.0 0.2 0.2 1.0, txt)]
        (return ())) SDLK_RETURN

checkArrived :: String -> String -> String -> Int -> [String] -> StateT SpaceState IO ()
checkArrived tgt planetname alleg pts str =
  when (tgt == planetname) (onActualArrival alleg pts str)


text:
{-# LANGUAGE OverloadedStrings #-}
module Data.Outputable.Instances where

import           Control.Applicative            (ZipList)
import           Data.Char                      (isAlphaNum)
import           Data.Complex                   (Complex)
import           Data.Functor.Identity          (Identity)
import           Data.Int                       (Int16, Int32, Int64, Int8)
import           Data.List.NonEmpty             (NonEmpty)
import qualified Data.Map                       as Map
import           Data.Monoid                    (All, Any, Dual, First, Last,
                                                 Product, Sum)
import           Data.Outputable.Class
import qualified Data.Semigroup                 as S (First, Last, Max, Min,
                                                      Option, WrappedMonoid)
import           Data.Text                      (Text, unpack)
import qualified Data.Text.Lazy                 as TL (Text, unpack)
import           Data.Version                   (Version)
import           Data.Void                      (Void)
import           Data.Word                      (Word16, Word32, Word64, Word8)
import           GHC.Generics
import           Numeric.Natural                (Natural)
import           System.Exit                    (ExitCode)
import           Text.PrettyPrint.HughesPJClass

viaPretty :: Pretty a => Int -> a -> Doc
viaPretty _ = pPrint

viaPrettyNum :: (Ord a, Num a, Pretty a) => Int -> a -> Doc
viaPrettyNum n x
  | n /= 0 && x < 0 = parens $ pPrint x
  | otherwise = pPrint x

viaShow _ x = text $ show x
viaShowNum n x
  | n /= 0 && x < 0 = parens $ viaShow n x
  | otherwise = viaShow n x

instance Outputable Double where pprPrec = viaPrettyNum
instance Outputable Float where pprPrec = viaPrettyNum
instance Outputable Int where pprPrec = viaPrettyNum
instance Outputable Int8 where pprPrec = viaShowNum
instance Outputable Int16 where pprPrec = viaShowNum
instance Outputable Int32 where pprPrec = viaShowNum
instance Outputable Int64 where pprPrec = viaShowNum
instance Outputable Integer where pprPrec = viaPrettyNum
instance Outputable Natural where pprPrec = viaShow
instance Outputable Word where pprPrec = viaShow
instance Outputable Word8 where pprPrec = viaShow
instance Outputable Word16 where pprPrec = viaShow
instance Outputable Word32 where pprPrec = viaShow
instance Outputable Word64 where pprPrec = viaShow
instance Outputable TL.Text where pprPrec n x = doubleQuotes $ text $ TL.unpack x
instance Outputable Text where pprPrec n x = doubleQuotes $ text $ unpack x
instance Outputable a => Outputable [a] where pprPrec _ = pprList

instance Outputable Char where
  pprPrec _ x = quotes $ char x
  pprList = pPrint

instance Outputable Bool
instance Outputable Ordering
instance Outputable ()
instance Outputable DecidedStrictness
instance Outputable SourceStrictness
instance Outputable SourceUnpackedness
instance Outputable Associativity
instance Outputable Fixity
instance Outputable Any
instance Outputable All
instance Outputable ExitCode
instance Outputable Version
instance Outputable Void
instance Outputable a => Outputable (Maybe a)
instance Outputable p => Outputable (Par1 p)
instance Outputable a => Outputable (NonEmpty a)
instance Outputable a => Outputable (Product a)
instance Outputable a => Outputable (Sum a)
instance Outputable a => Outputable (Dual a)
instance Outputable a => Outputable (Last a)
instance Outputable a => Outputable (First a)
instance Outputable a => Outputable (Identity a)
instance Outputable a => Outputable (ZipList a)
instance Outputable a => Outputable (S.Option a)
instance Outputable m => Outputable (S.WrappedMonoid m)
instance Outputable a => Outputable (S.Last a)
instance Outputable a => Outputable (S.First a)
instance Outputable a => Outputable (S.Max a)
instance Outputable a => Outputable (S.Min a)
instance Outputable a => Outputable (Complex a)
instance (Outputable a, Outputable b) => Outputable (Either a b)
instance (Outputable a, Outputable b) => Outputable (a, b) where
  pprPrec _ (a, b) = parens $ fsep $ punctuate comma [ppr a, ppr b]
instance (Outputable a, Outputable b, Outputable c) => Outputable (a, b, c) where
  pprPrec _ (a, b, c) = parens $ fsep $ punctuate comma [ppr a, ppr b, ppr c]
instance (Outputable a, Outputable b, Outputable c, Outputable d) => Outputable (a, b, c, d) where
  pprPrec _ (a, b, c, d) = parens $ fsep $ punctuate comma [ppr a, ppr b, ppr c, ppr d]
instance (Outputable a, Outputable b, Outputable c, Outputable d, Outputable e) => Outputable (a, b, c, d, e) where
  pprPrec _ (a, b, c, d, e) = parens $ fsep $ punctuate comma [ppr a, ppr b, ppr c, ppr d, ppr e]

instance (Outputable k, Outputable v) => Outputable (Map.Map k v) where
  pprPrec _ m = "fromList" <+> ppr (Map.toList m)


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE DefaultSignatures, TypeOperators, FlexibleContexts, TypeSynonymInstances, FlexibleInstances #-}
module Pure.Data.Default where

import Control.Applicative

import Data.Complex
import Data.Int
import Data.Monoid
import Data.Ratio
import Data.Word

import GHC.Generics

class Default a where
  def :: a
  default def :: (Generic a, GDefault (Rep a)) => a
  def = to gdef

instance Default Bool where def = False -- like other languages with construction defaults

instance Default () where def = ()
instance Default Ordering where def = EQ
instance Default Any where def = Any False
instance Default All where def = All True
instance Default (Last a) where def = Last Nothing
instance (Num a) => Default (Sum a) where def = Sum 0
instance (Num a) => Default (Product a) where def = Product 1
instance Default (Endo a) where def = Endo id
instance Default a => Default (Const a b) where def = Const def
instance Default (Maybe a) where def = Nothing
instance Default [a] where def = []

-- Note that (def :: Try a) /= mempty
-- instance Default (Try a) where def = Trying

-- instance Default Txt where def = mempty

-- instance (Eq a, Hashable a) => Default (HashMap.HashMap a b) where
--   def = mempty

-- instance Default Value where
--   def =
-- #ifdef __GHCJS__
--     nullValue
-- #else
--     Null
-- #endif

-- instance Default Obj where def = mempty

-- instance Default Micros where def = 0
-- instance Default Millis where def = 0
instance Default Int where def = 0
instance Default Int8 where def = 0
instance Default Int16 where def = 0
instance Default Int32 where def = 0
instance Default Int64 where def = 0
instance Default Word where def = 0
instance Default Word8 where def = 0
instance Default Word16 where def = 0
instance Default Word32 where def = 0
instance Default Word64 where def = 0
instance Default Integer where def = 0
instance Default Float where def = 0
instance Default Double where def = 0
instance (Integral a) => Default (Ratio a) where def = 0
instance (Default a,RealFloat a) => Default (Complex a) where def = def :+ def
instance {-# OVERLAPPABLE #-} Default r => Default (x -> r) where def = const def
instance {-# OVERLAPPING #-} Default (a -> a) where def = id
instance Default a => Default (IO a) where def = return def
instance Default a => Default (Dual a) where def = Dual def

instance (Default a, Default b) => Default (a, b) where def = (def, def)
instance (Default a, Default b, Default c) => Default (a, b, c) where def = (def, def, def)
instance (Default a, Default b, Default c, Default d) => Default (a, b, c, d) where def = (def, def, def, def)
instance (Default a, Default b, Default c, Default d, Default e) => Default (a, b, c, d, e) where def = (def, def, def, def, def)
instance (Default a, Default b, Default c, Default d, Default e, Default f) => Default (a, b, c, d, e, f) where def = (def, def, def, def, def, def)
instance (Default a, Default b, Default c, Default d, Default e, Default f, Default g) => Default (a, b, c, d, e, f, g) where def = (def, def, def, def, def, def, def)

-- Inspired by Lukas Mai's data-default-class with a default instance for
-- sum types based on lexicographical order - similar to the Ord and Enum
-- instances
class GDefault f where
  gdef :: f a

instance GDefault V1 where
  gdef = undefined

instance GDefault U1 where
  gdef = U1

instance (Default a) => GDefault (K1 i a) where
  gdef = K1 def

instance (GDefault a, GDefault b) => GDefault (a :*: b) where
  gdef = gdef :*: gdef

instance (GDefault a, GDefault b) => GDefault (a :+: b) where
  gdef = L1 gdef

instance (GDefault a) => GDefault (M1 i c a) where
  gdef = M1 gdef

max_stars_repo_name: grumply/pure-default

max_issues_repo_name: grumply/pure-default

max_forks_repo_name: grumply/pure-default


text:
module Main where

import Test.Hspec
import Data.Word
import Data.Complex
import Data.Vector.Storable (Vector)
import Data.Vector.Storable as SV
import Codec.Picture
import Data.Monoid
import Control.Monad.State
import Numeric.LinearAlgebra hiding ((<>))
import Test.Hspec.QuickCheck
import Test.QuickCheck.Property
import Test.QuickCheck

import FFT
import ArrayPoking
import RealMain
import Algoritmo

main :: IO ()
main = hspec $ do
  describe "FFT" $ do
    it "Valor conocido de FFT" $
      fft ex_vec1 `shouldBe` SV.fromList [2, 4 :+ 4, -6, 4 :+ (-4)]

    it "ifft . fft = id" $
      ifft (fft ex_vec1) `shouldBe` ex_vec1

  describe "array poking" $ do
    it "Separa inicio" $ do
      lena_img <- loadImage lena_path
      let red   = extractColor 0 lena_img
          green = extractColor 1 lena_img
          blue  = extractColor 2 lena_img
      SV.take 3 (imageData lena_img) `shouldBe`
        (SV.take 1 red <> SV.take 1 green <> SV.take 1 blue)

    it "Recombinacion" $ do
      lena_img <- loadImage lena_path
      let red   = extractColor 0 lena_img
          green = extractColor 1 lena_img
          blue  = extractColor 2 lena_img
      imageData lena_img `shouldBe` combineColors red green blue

    it "Recombinando %error" $ do
      -- lena_img@(Image width height _) <- loadImage lena_path
      let alt = SV.fromList [4, 2, 3, 7, 10, 8, 10, 14]
          res = evalState (secAnalitica 2 alt >> secSintetica 2) []
      norm_2 (res - alt) / norm_2 alt < 0.1 `shouldBe` True

  -- prop "Recombinando %error" $
  --     -- lena_img@(Image width height _) <- loadImage lena_path
  --     let res alt = evalState (secAnalitica 4 alt >> secSintetica) []
  --         -- alt = SV.fromList [4, 2, 3, 7, 10, 8, 10, 14, 40, 22, 90, 33, 11, 34, 38, 50]
  --     in  (forAll (fmap SV.fromList arbitrary)
  --         (\a -> norm_2 (res a - a) / norm_2 a < 0.01))

  -- describe "Pruebas en Lena"
  --   it "Recombinacion Lena" $ do
  --     lena_img <- loadImage lena_path
  --     let red   = extractColor 0 lena_img
  --         green = extractColor 1 lena_img
  --         blue  = extractColor 2 lena_img
  --         go color = evalState (secAnalitica color >>= secSintetica)
  --                              (VS.length color)
  --     imageData lena_img `shouldBe` combineColors red green blue

{- Ejemplo libro pagina 126 -}
ex_vec1 :: Vector (Complex Double)
ex_vec1 = SV.fromList [ 1, 0, -3, 4]

max_stars_repo_name: RubenAstudillo/proyecto-lena

max_issues_repo_name: RubenAstudillo/proyecto-lena

max_forks_repo_name: RubenAstudillo/proyecto-lena


text:
module Main (main) where

import AI.HNN.FF.Network
import Data.List (maximumBy)
import Data.Ord (comparing)
import Data.Word (Word8)
import MnistLoader
import Numeric.LinearAlgebra
import Text.Printf

main :: IO ()
main = do
  (samples, validation) <- load_data_wrapper "../data/"
  net <- createNetwork 784 [30] 10
  let net' = trainNTimes 30 3.0 sigmoid sigmoid' net samples
      result = map (\(a, b) -> (output net' sigmoid a, b)) validation
  printf "%.2f%% correct\n" (100 * validate result)

validate :: [(Vector Double, Word8)] -> Float
validate results =
  correct / (fromIntegral $ length results)
  where
    correct = sum [ 1 | (x, y) <- results, maxListIndex x == y]
    maxListIndex :: Vector Double -> Word8
    maxListIndex = fst . maxPair . toList
    maxPair :: [Double] -> (Word8, Double)
    maxPair = maximumBy (comparing snd) . zip [0..]

max_stars_repo_name: KaiHa/NNDigitReader

max_issues_repo_name: KaiHa/NNDigitReader

max_forks_repo_name: KaiHa/NNDigitReader


text:
{-# LANGUAGE TypeApplications #-}

module Main where

import Control.Monad.State
import Control.Monad.Trans.Class
import Numeric.LinearAlgebra.HMatrix hiding (scale)
import qualified Graphics.Rendering.OpenGL as GL
import Foreign.Storable
import GLCode
import Sketch
import Points
import Lines

v2f :: Real a => a -> a -> GL.Vertex2 Float
v2f x y = GL.Vertex2 (realToFrac x) (realToFrac y)

translate :: (Fractional a, Storable a) => a -> a -> a -> Matrix a
translate x y z = (4 >< 4) [ 1.0, 0.0, 0.0, 0.0,
                             0.0, 1.0, 0.0, 0.0,
                             0.0, 0.0, 1.0, 0.0,
                             x  , y  , z  , 1.0 ]

scale :: (Fractional a, Storable a) => a -> a -> a -> Matrix a
scale x y z = (4 >< 4) [ x  , 0.0, 0.0, 0.0,
                         0.0, y  , 0.0, 0.0,
                         0.0, 0.0, z  , 0.0,
                         0.0, 0.0, 0.0, 1.0 ]

rotate :: (Fractional a, Storable a, Floating a) => a -> Matrix a
rotate theta = let c = cos theta
                   s = sin theta
               in (4 >< 4) [ c, -s, 0.0, 0.0,
                             s, c, 0.0, 0.0,
                             0.0, 0.0, 1.0, 0.0,
                             0.0, 0.0, 0.0, 1.0]

translateBy :: Float -> Float -> Float -> StateT (Matrix Float) (StateT World IO) ()
translateBy x y z = modify (translate x y z `mul`)

scaleBy :: Float -> StateT (Matrix Float) (StateT World IO) ()
scaleBy s = modify (scale s s s `mul`)

rotateBy :: Float -> StateT (Matrix Float) (StateT World IO) ()
rotateBy s = modify (rotate s `mul`)

main :: IO ()
main = mainLoopState 512 512 1 (ident @Float 4) $ \time -> do

        put (ident @Float 4)

        let t = 5*realToFrac time :: Float

        GL.clearColor GL.$= GL.Color4 0.0 0.0 0.0 1
        io $ GL.clear [GL.ColorBuffer]

        GL.lineSmooth GL.$= GL.Enabled
        GL.lineWidth GL.$= 2
        GL.hint GL.LineSmooth GL.$= GL.Nicest

        modify (rotate (0.1*time) `mul`)
        replicateM_ 23 $ do
            let r = 10+5*sin (0.71*time)
            rotateBy (pi/r)
            let s = sin (0.9*time)
            scaleBy s
            translateBy 0.01 0.01 0.0
            transform <- get

            lift $ setUniform "turtle"
                       "transform" transform
            lift $ drawLine "turtle"
                     "vPosition" [v2f (-1.0) 0.0, v2f 1.0 0.0]

            io GL.flush

max_stars_repo_name: dpiponi/SketchBox

max_issues_repo_name: dpiponi/SketchBox

max_forks_repo_name: dpiponi/SketchBox


text:
module MachineLearning.LogisticModelTest
(
  tests
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Approx
import Test.HUnit.Plus

import MachineLearning.DataSets (dataset2)

import qualified Numeric.LinearAlgebra as LA
import qualified MachineLearning as ML
import MachineLearning.Optimization (checkGradient)
import MachineLearning.Model
import MachineLearning.LogisticModel
import MachineLearning.Regularization (Regularization(..))

(x, y) = ML.splitToXY dataset2

x1 = ML.addBiasDimension $ ML.mapFeatures 6 x
onesTheta :: LA.Vector LA.R
onesTheta = LA.konst 1 (LA.cols x1)
zeroTheta :: LA.Vector LA.R
zeroTheta = LA.konst 0 (LA.cols x1)

gradientCheckingEps = 1e-3

checkGradientTest lambda theta = do
  let diffs = take 5 $ map (\e -> checkGradient Logistic lambda x1 y theta e) [1e-3, 1.1e-3 ..]
      diff = minimum diffs
  assertApproxEqual (show theta) gradientCheckingEps 0 diff

tests = [ testGroup "sigmoid" [
            testCase "zero" $ assertApproxEqual "" 1e-10 0.5 (sigmoid 0)
            , testCase "big positive value" $ assertApproxEqual "" 1e-10 1 (sigmoid 10e10)
            , testCase "big negative value" $ assertApproxEqual "" 1e-10 0 (sigmoid $ -10e10)
            , testCase "quite big positive value" $ assertApproxEqual "" 1e-2 1 (sigmoid 100)
            , testCase "quite big negative value" $ assertApproxEqual "" 1e-2 0 (sigmoid $ -100)
            ]
          , testGroup "sigmoidGradient" [
              testCase "zero" $ assertApproxEqual "" 1e-10 0.25 (sigmoidGradient 0)
              , testCase "big positive value" $ assertApproxEqual "" 1e-10 0 (sigmoidGradient 10e10)
              , testCase "big negative value" $ assertApproxEqual "" 1e-10 0 (sigmoidGradient $ -10e10)
              , testCase "quite big positive value" $ assertApproxEqual "" 1e-2 0 (sigmoidGradient 100)
              , testCase "quite big negative value" $ assertApproxEqual "" 1e-2 0 (sigmoidGradient $ -100)
              , testCase "small positive value" $ assertApproxEqual "" 1e-2 0.2 (sigmoidGradient 1)
              , testCase "small negative value" $ assertApproxEqual "" 1e-2 0.2 (sigmoidGradient $ -1)
              ]
          , testGroup "model" [
              testCase "cost, lambda = 0" $ assertApproxEqual "" 1e-3 2.020 (cost Logistic (L2 0) x1 y onesTheta)
              , testCase "cost, lambda = 1" $ assertApproxEqual "" 1e-3 2.135 (cost Logistic (L2 1) x1 y onesTheta)
              , testCase "cost, lambda = 1000" $ assertApproxEqual "" 1e-3 116.427 (cost Logistic (L2 1000) x1 y onesTheta)
              , testCase "gradient, lambda = 0" $ assertVector "" 1e-5 gradient_l0 (gradient Logistic (L2 0) x1 y onesTheta)
              , testCase "gradient, lambda = 1" $ assertVector "" 1e-5 gradient_l1 (gradient Logistic (L2 1) x1 y onesTheta)
              , testCase "gradient, lambda = 1000" $ assertVector "" 1e-5 gradient_l1000 (gradient Logistic (L2 1000) x1 y onesTheta)
              ]
          , testGroup "gradient checking" [
              testCase "non-zero theta, non-zero lambda" $ checkGradientTest (L2 2) onesTheta
              , testCase "zero theta, non-zero lambda" $ checkGradientTest (L2 2) zeroTheta
              , testCase "non-zero theta, zero lambda" $ checkGradientTest (L2 0) onesTheta
              , testCase "zero theta, zero lambda" $ checkGradientTest (L2 0) zeroTheta
              , testCase "non-zero theta, no reg" $ checkGradientTest RegNone onesTheta
              , testCase "zero theta, no reg" $ checkGradientTest RegNone zeroTheta

              ]
        ]

gradient_l0 = LA.vector [0.34604507367924525,7.660615656904722e-2,0.11004999290013262,0.14211701951318526,7.4399123914273965e-3,0.15963981400206023,5.864636026438637e-2,2.369595228045015e-2,1.7568631688383615e-2,9.87226947583087e-2,8.878427090266403e-2,2.509755728155993e-3,3.348199373717313e-2,1.0975419586624715e-3,0.11520318246520422,5.0480764463147594e-2,1.0229511332420786e-2,8.818652179760128e-3,1.5052078581608905e-2,6.655810974747264e-3,9.010665405440292e-2,6.480865498500844e-2,2.039892642614106e-3,1.4231095091583643e-2,5.737477462013599e-4,1.71609000731311e-2,-2.437841009425525e-4,9.753746346629336e-2]

gradient_l1 = LA.vector [0.34604507367924525,8.508073284023365e-2,0.11852456917131905,0.15059159578437167,1.5914488662613836e-2,0.16811439027324665,6.71209365355728e-2,3.217052855163659e-2,2.6043207959570054e-2,0.10719727102949513,9.725884717385046e-2,1.0984331999342433e-2,4.1956570008359576e-2,9.572118229848912e-3,0.12367775873639067,5.895534073433403e-2,1.8704087603607224e-2,1.729322845094657e-2,2.3526654852795346e-2,1.5130387245933704e-2,9.858123032558937e-2,7.328323125619488e-2,1.0514468913800546e-2,2.270567136277008e-2,9.048324017387801e-3,2.563547634431754e-2,8.230792170243889e-3,0.10601203973747979]

gradient_l1000 = LA.vector [0.34604507367924525,8.551182427755489,8.584626264086573,8.616693290699626,8.48201618357787,8.6342160851885,8.533222631450828,8.498272223466891,8.492144902874823,8.57329896594475,8.563360542089104,8.477086026914597,8.508058264923614,8.475673813145104,8.589779453651644,8.525057035649588,8.484805782518862,8.483394923366202,8.489628349768049,8.481232082161188,8.564682925240843,8.539384926171449,8.476616163829055,8.488807366278024,8.475150018932643,8.491737171259572,8.474332487085498,8.572113734652733]

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ExplicitNamespaces #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Presburger #-}

{- | Algorithms for zero-dimensional ideals.

   Since 0.4.0.0
-}
module Algebra.Algorithms.ZeroDim
  ( -- * Root finding for zero-dimensional ideal
    solveM,
    solve',
    solveViaCompanion,
    solveLinear,

    -- * Radical computation
    radical,
    isRadical,

    -- * Converting monomial ordering to Lex using FGLM algorithm
    fglm,
    fglmMap,

    -- ** Internal helper function
    solveWith,
    univPoly,
    reduction,
    matrixRep,
    subspMatrix,
    vectorRep,
  )
where

import Algebra.Algorithms.FGLM
import Algebra.Algorithms.Groebner
import Algebra.Instances ()
import Algebra.Internal hiding (OLt)
import qualified Algebra.Matrix as AM
import Algebra.Prelude.Core
import Algebra.Ring.Polynomial.Quotient
import Control.Lens
import Control.Monad.Loops (whileM_)
import Control.Monad.Random hiding (next)
import Control.Monad.Reader (runReaderT)
import Control.Monad.ST (runST)
import Data.Complex (Complex (..), magnitude)
import Data.Convertible (Convertible, convert)
import qualified Data.Matrix as M
import Data.Maybe (fromJust)
import Data.Reflection (Reifies)
import Data.STRef.Strict (newSTRef)
import qualified Data.Sized as SV
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV
import qualified Numeric.Algebra as NA
import qualified Numeric.LinearAlgebra as LA
import qualified Prelude as P

{- | Finds complex approximate  roots of given zero-dimensional ideal,
   using randomized altorithm.

   See also @'solve''@ and @'solveViaCompanion'@.
-}
solveM ::
  forall m r ord n.
  ( Normed r
  , Ord r
  , MonadRandom m
  , Field r
  , CoeffRing r
  , KnownNat n
  , IsMonomialOrder n ord
  , Convertible r Double
  , 0 < n
  ) =>
  Ideal (OrderedPolynomial r ord n) ->
  m [Sized n (Complex Double)]
solveM ideal =
  {-# SCC "solveM" #-}
  withKnownNat (sSucc (sNat :: SNat n)) $
    reifyQuotient (radical ideal) $ \pxy ->
      case standardMonomials' pxy of
        Just bs -> step 10 (length bs)
        Nothing -> error "Given ideal is not zero-dimensional"
  where
    step bd len =
      {-# SCC "solveM/step" #-}
      do
        coeffs <-
          {-# SCC "solveM/coeff-gen" #-}
          replicateM (sNatToInt (sSucc (sNat :: SNat n))) $ getRandomR (- bd, bd)
        let vars = one : SV.toList allVars
            f = sum $ zipWith (.*.) (map (NA.fromInteger :: Integer -> r) coeffs) vars
        case solveWith f ideal of
          Nothing -> step (bd * 2) len
          Just sols -> return sols

{- | @solveWith f is@ finds complex approximate roots of the given zero-dimensional @n@-variate polynomial system @is@,
   using the given relatively prime polynomial @f@.
-}
solveWith ::
  forall r n ord.
  ( DecidableZero r
  , Normed r
  , Ord r
  , Field r
  , CoeffRing r
  , 0 < n
  , IsMonomialOrder n ord
  , KnownNat n
  , Convertible r Double
  ) =>
  OrderedPolynomial r ord n ->
  Ideal (OrderedPolynomial r ord n) ->
  Maybe [Sized n (Complex Double)]
solveWith f0 i0 =
  {-# SCC "solveWith" #-}
  withKnownNat (sSucc (sNat :: SNat n)) $
    reifyQuotient (radical i0) $ \pxy ->
      let ideal = gBasis' pxy
          Just base = map (leadingMonomial . quotRepr) <$> standardMonomials' pxy
       in case {-# SCC "findOne" #-} elemIndex one base of
            Nothing -> Just []
            Just cind ->
              let f = modIdeal' pxy f0
                  vars =
                    sortBy (flip $ comparing snd) $
                      map (\on -> (on, leadingMonomial $ var on `asTypeOf` f0)) $
                        enumOrdinal $ sArity' f0
                  inds = flip map vars $
                    second $ \b ->
                      case findIndex (== b) base of
                        Just ind -> Right ind
                        Nothing ->
                          let Just g = find ((== b) . leadingMonomial) ideal
                              r = leadingCoeff g
                              answer = mapCoeff toComplex $ injectCoeff (recip r) * (toPolynomial (leadingTerm g) - g)
                           in Left answer
                  mf = AM.fromLists $ map (map toComplex) $ matrixRep f
                  (_, evecs) = LA.eig $ LA.tr mf
                  calc vec =
                    {-# SCC "calc" #-}
                    let c = vec LA.! cind
                        phi (idx, Right nth) acc = acc & ix idx .~ (vec LA.! nth) P./ c
                        phi (idx, Left g) acc = acc & ix idx .~ substWith (*) acc g
                     in if c == 0
                          then Nothing
                          else Just $ foldr ({-# SCC "rewrite-answer" #-} phi) (SV.replicate (sArity' f0) (error "indec!")) inds
               in sequence $ map calc $ LA.toColumns evecs

{- | @'solve'' err is@ finds numeric approximate root of the
   given zero-dimensional polynomial system @is@,
   with error <@err@.

   See also @'solveViaCompanion'@ and @'solveM'@.
-}
solve' ::
  forall r n ord.
  ( Field r
  , CoeffRing r
  , KnownNat n
  , 0 < n
  , IsMonomialOrder n ord
  , Convertible r Double
  ) =>
  Double ->
  Ideal (OrderedPolynomial r ord n) ->
  [Sized n (Complex Double)]
solve' err ideal =
  reifyQuotient ideal $ \ii ->
    if gBasis' ii == [one]
      then []
      else
        let vs =
              map (nub . LA.toList . LA.eigenvalues . AM.fromLists . map (map toComplex) . matrixRep . modIdeal' ii) $
                SV.toList allVars
            mul p q = toComplex p * q
         in [ xs
            | xs0 <- sequence vs
            , let xs = SV.unsafeFromList' xs0
            , all ((< err) . magnitude . substWith mul xs) $ generators ideal
            ]

{- | Given a zero-dimensional ideal \(I\),
   \('subspMatrix' i I\) computes a multiplication matrix
   \(m_{x_i} \mathbin{\upharpoonright} V_i \) by \(x_i\)
   restricted to the subspace \(V_i = \mathop{\mathrm{span}}(\left\{\ x_i^n \ \middle|\ 1 \leq n \right\})\) of \(k[\mathbf{X}]\).
-}
subspMatrix ::
  forall r n ord.
  (Show r, Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord) =>
  Ordinal n ->
  Ideal (OrderedPolynomial r ord n) ->
  M.Matrix r
subspMatrix on ideal =
  let poly = univPoly on ideal
      v = var on :: OrderedPolynomial r ord n
      dim = totalDegree' poly
      cfs = [negate $ coeff (leadingMonomial $ pow v (j :: Natural)) poly | j <- [0 .. fromIntegral (dim - 1)]]
      leftTops = M.fromLists [replicate (dim - 1) zero]
      leftBots = fmap unwrapAlgebra (M.identity (dim - 1))
      lfts =
        leftTops
          M.<-> leftBots
      rights = M.colVector (V.fromList cfs)
   in lfts M.<|> rights

{- | @'solveViaCompanion' err is@ finds numeric approximate root of the
   given zero-dimensional polynomial system @is@,
   with error <@err@.

   See also @'solve''@ and @'solveM'@.
-}
solveViaCompanion ::
  forall r ord n.
  (Show r, Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord, Convertible r Double) =>
  Double ->
  Ideal (OrderedPolynomial r ord n) ->
  [Sized n (Complex Double)]
solveViaCompanion err ideal =
  if calcGroebnerBasis ideal == [one]
    then []
    else
      let vs =
            map (nub . LA.toList . LA.eigenvalues . LA.fromLists . matToLists . fmap toComplex . flip subspMatrix ideal) $
              enumOrdinal (sNat :: SNat n)
          mul p q = toComplex p * q
       in [ xs
          | xs0 <- sequence vs
          , let xs = SV.unsafeFromList' xs0
          , all ((< err) . magnitude . substWith mul xs) $ generators ideal
          ]

matToLists :: M.Matrix a -> [[a]]
matToLists mat = [V.toList $ M.getRow i mat | i <- [1 .. M.nrows mat]]

matrixRep ::
  ( DecidableZero t
  , Eq t
  , Field t
  , KnownNat n
  , IsMonomialOrder n order
  , Reifies ideal (QIdeal (OrderedPolynomial t order n))
  ) =>
  Quotient (OrderedPolynomial t order n) ideal ->
  [[t]]
matrixRep f =
  {-# SCC "matrixRep" #-}
  case standardMonomials of
    Just [] -> []
    Just bases ->
      let anss = map (quotRepr . (f *)) bases
       in transpose $ map (\a -> map (flip coeff a . leadingMonomial . quotRepr) bases) anss
    Nothing -> error "Not finite dimension"

toComplex :: Convertible a Double => a -> Complex Double
toComplex a = convert a :+ 0

-- | Calculates n-th reduction of f: @f `div` < f, ∂_{x_n} f >@.
reduction ::
  (CoeffRing r, KnownNat n, IsMonomialOrder n ord, Field r) =>
  Ordinal n ->
  OrderedPolynomial r ord n ->
  OrderedPolynomial r ord n
reduction on f =
  {-# SCC "reduction" #-}
  let df = {-# SCC "differentiate" #-} diff on f
   in snd $ head $ f `divPolynomial` calcGroebnerBasis (toIdeal [f, df])

-- | Calculate the monic generator of \( k[X_0, ..., X_n] \cap k[X_i]\).
univPoly ::
  forall r ord n.
  (Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord) =>
  Ordinal n ->
  Ideal (OrderedPolynomial r ord n) ->
  OrderedPolynomial r ord n
univPoly nth ideal =
  {-# SCC "univPoly" #-}
  reifyQuotient ideal $ \pxy ->
    if gBasis' pxy == [one]
      then one
      else
        let x = var nth
            monomDeg =
              fromIntegral $
                length $
                  fromJust $
                    standardMonomials' pxy
            p0 : pows =
              [ fmap WrapAlgebra $ vectorRep $ modIdeal' pxy (pow x i)
              | i <- [0 :: Natural .. monomDeg]
              ]
            step m ~(p : ps) =
              {-# SCC "univPoly/step" #-}
              case solveLinear m p of
                Nothing -> {-# SCC "recur" #-} step ({-# SCC "consCol" #-} m M.<|> M.colVector p) ps
                Just ans ->
                  let cur = fromIntegral $ V.length ans :: Natural
                   in {-# SCC "buildRelation" #-}
                      pow x cur
                        - sum
                          ( zipWith
                              (.*.)
                              (fmap unwrapAlgebra $ V.toList ans)
                              [pow x i | i <- [0 :: Natural .. cur P.- 1]]
                          )
         in step (M.colVector p0) pows

-- | Solves linear system. If the given matrix is degenerate, this returns @Nothing@.
solveLinear ::
  (Ord r, P.Fractional r) =>
  M.Matrix r ->
  V.Vector r ->
  Maybe (V.Vector r)
solveLinear mat vec =
  {-# SCC "solveLinear" #-}
  if ({-# SCC "uRank" #-} uRank u) < uRank u' || M.diagProd u == 0 || uRank u < M.ncols mat
    then Nothing
    else
      let ans = M.getCol 1 $ p P.* M.colVector vec
          lsol = {-# SCC "solveL" #-} solveL ans
          cfs = M.getCol 1 $ q P.* M.colVector ({-# SCC "solveU" #-} solveU lsol)
       in Just cfs
  where
    Just (u, l, p, q, _, _) = M.luDecomp' mat
    Just (u', _, _, _, _, _) = M.luDecomp' (mat M.<|> M.colVector vec)
    uRank = V.foldr (\a acc -> if a /= 0 then acc + 1 else acc) (0 :: Int) . M.getDiag
    solveL v = V.create $ do
      let stop = min (M.ncols l) (M.nrows l)
      mv <- MV.replicate (M.ncols l) 0
      forM_ [0 .. stop - 1] $ \i -> do
        MV.write mv i $ v V.! i
        forM_ [0, 1 .. min (i -1) (M.ncols l - 1)] $ \j -> do
          a <- MV.read mv i
          b <- MV.read mv j
          MV.write mv i $ a P.- (l M.! (i + 1, j + 1)) P.* b
      return mv
    solveU v = V.create $ do
      let stop = min (M.ncols u) (M.nrows u)
      mv <- MV.replicate (M.ncols u) 0
      forM_ [stop - 1, stop - 2 .. 0] $ \i -> do
        MV.write mv i $ v V.! i
        forM_ [i + 1, i + 2 .. M.ncols u -1] $ \j -> do
          a <- MV.read mv i
          b <- MV.read mv j
          MV.write mv i $ a P.- (u M.! (i + 1, j + 1)) P.* b
        a0 <- MV.read mv i
        MV.write mv i $ a0 P./ (u M.! (i + 1, i + 1))
      return mv

-- | Calculate the radical of the given zero-dimensional ideal.
radical ::
  forall r ord n.
  ( Ord r
  , CoeffRing r
  , KnownNat n
  , Field r
  , IsMonomialOrder n ord
  ) =>
  Ideal (OrderedPolynomial r ord n) ->
  Ideal (OrderedPolynomial r ord n)
radical ideal =
  {-# SCC "radical" #-}
  let gens = {-# SCC "calcGens" #-} map (\on -> reduction on $ univPoly on ideal) $ enumOrdinal (sNat :: SNat n)
   in toIdeal $ calcGroebnerBasis $ toIdeal $ generators ideal ++ gens

-- | Test if the given zero-dimensional ideal is radical or not.
isRadical ::
  forall r ord n.
  ( Ord r
  , CoeffRing r
  , KnownNat n
  , 0 < n
  , Field r
  , IsMonomialOrder n ord
  ) =>
  Ideal (OrderedPolynomial r ord n) ->
  Bool
isRadical ideal =
  let gens =
        map (\on -> reduction on $ univPoly on ideal) $
          enumOrdinal (sNat :: SNat n)
   in all (`isIdealMember` ideal) gens

-- * FGLM

{- | Calculate the Groebner basis w.r.t. lex ordering of the zero-dimensional ideal using FGLM algorithm.
   If the given ideal is not zero-dimensional this function may diverge.
-}
fglm ::
  ( Ord r
  , KnownNat n
  , Field r
  , IsMonomialOrder n ord
  , 0 < n
  ) =>
  Ideal (OrderedPolynomial r ord n) ->
  ([OrderedPolynomial r Lex n], [OrderedPolynomial r Lex n])
fglm ideal = reifyQuotient ideal $ \pxy ->
  fglmMap (vectorRep . modIdeal' pxy)

-- | Compute the kernel and image of the given linear map using generalized FGLM algorithm.
fglmMap ::
  forall k ord n.
  ( Ord k
  , Field k
  , 0 < n
  , IsMonomialOrder n ord
  , CoeffRing k
  , KnownNat n
  ) =>
  -- | Linear map from polynomial ring.
  (OrderedPolynomial k ord n -> V.Vector k) ->
  -- | The tuple of:
  --
  --     * lex-Groebner basis of the kernel of the given linear map.
  --
  --     * The vector basis of the image of the linear map.
  ( [OrderedPolynomial k Lex n]
  , [OrderedPolynomial k Lex n]
  )
fglmMap l = runST $ do
  env <- FGLMEnv l <$> newSTRef [] <*> newSTRef [] <*> newSTRef Nothing <*> newSTRef one
  flip runReaderT env $ do
    mainLoop
    whileM_ toContinue $ nextMonomial >> mainLoop
    (,) <$> look gLex <*> (map (changeOrder Lex) <$> look bLex)

mainLoop ::
  (DecidableZero r, Ord r, KnownNat n, Field r, IsMonomialOrder n o) =>
  Machine s r o n ()
mainLoop = do
  m <- look monomial
  let f = toPolynomial (one, changeMonomialOrderProxy Proxy m)
  lx <- image f
  bs <- mapM image =<< look bLex
  let mat = foldr (M.<|>) (M.fromList 0 0 []) $ map (M.colVector . fmap WrapAlgebra) bs
      cond
        | null bs =
          if V.all (== zero) lx
            then Just $ V.replicate (length bs) 0
            else Nothing
        | otherwise = solveLinear mat (fmap WrapAlgebra lx)
  case cond of
    Nothing -> do
      proced .== Nothing
      bLex @== (f :)
    Just cs -> do
      bps <- look bLex
      let g = changeOrder Lex $ f - sum (zipWith (.*.) (V.toList $ fmap unwrapAlgebra cs) bps)
      proced .== Just (changeOrder Lex f)
      gLex @== (g :)

toContinue ::
  forall s r o n.
  ( 0 < n
  , Ord r
  , KnownNat n
  , Field r
  ) =>
  Machine s r o n Bool
toContinue = do
  mans <- look proced
  case mans of
    Nothing -> return True
    Just g -> do
      let xLast = P.maximum allVars `asTypeOf` g
      return $ not $ leadingMonomial g `isPowerOf` leadingMonomial xLast

nextMonomial ::
  forall s r ord n.
  (CoeffRing r, KnownNat n) =>
  Machine s r ord n ()
nextMonomial = do
  m <- look monomial
  gs <- map leadingMonomial <$> look gLex
  let next =
        fst $
          maximumBy
            (comparing snd)
            [ (OrderedMonomial monom, fromEnum od)
            | od <- enumOrdinal (sNat :: SNat n)
            , let monom = beta (getMonomial m) od
            , all (not . (`divs` OrderedMonomial monom)) gs
            ]
  monomial .== next

beta :: KnownNat n => Monomial n -> Ordinal n -> Monomial n
beta xs o@(OLt k) =
  let n = sizedLength xs
   in (SV.take (sSucc k) $ xs & ix o +~ 1) SV.++ SV.replicate (n %- sSucc k) 0
beta _ _ = error "beta: Bug in ghc!"

max_stars_repo_name: hangingman/computational-algebra

max_issues_repo_name: hangingman/computational-algebra

max_forks_repo_name: hangingman/computational-algebra


text:
#!/usr/bin/env runhaskell

import Data.Complex
import System.Environment

import Escape

usage :: IO ()
usage = do
    progName <- getProgName
    putStrLn $ "usage: " ++ progName ++ " <fractal type> <args>"
    putStrLn $ "  fractal type: mandelbrot"
    putStrLn $ "                julia x y (x y contructs a complex number x+yi)"
    putStrLn $ "                ship"
    putStrLn $ "                nova p r (f(z) = z^p - 1, z' = z - r*f(z)/f'(z))"

main = do
    (op:args) <- fmap (++ [""] ++ [""]) getArgs -- make sure it matches pattern (_:_)
    case op of
        "help"       -> usage
        "mandelbrot" -> mandelbrot
        "julia"      -> do
            let x = read $ args !! 0 :: Double
            let y = read $ args !! 1 :: Double
            julia (x :+ y)
        "ship"       -> burningShip
        "nova"       -> do
            let p = read $ args !! 0 :: Int
            let r = read $ args !! 1 :: Double
            nova p r
        _            -> usage


text:
{-# LANGUAGE CPP              #-}
{-# LANGUAGE DataKinds        #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs            #-}
{-# LANGUAGE TemplateHaskell  #-}
{-# LANGUAGE TypeOperators    #-}

module Test.Grenade.Sys.Networks where

import qualified Numeric.LinearAlgebra        as LA
import qualified Numeric.LinearAlgebra.Data   as D
import           Numeric.LinearAlgebra.Static (L)
import qualified Numeric.LinearAlgebra.Static as H

import           Hedgehog

import           Test.Grenade.Sys.Utils
import           Test.Hedgehog.Compat

import           Data.List
import           Data.Maybe
import           Data.Function
import           Data.Either
import           System.FilePath

import           Grenade
import           Grenade.Utils.PascalVoc

-- Checks that the differences in the reference and the output from the network are "acceptable". Due to floating
-- point errors, using doubles instead of floats, etc. the outputs can be not precisely the same. We check that
-- each box given by the reference output has a corresponding box in ours which roughly overlaps with it, detects
-- the same object, with approximately the same confidence, and vice versa for our output compared to the reference.
yoloAcceptableDifference :: [(Int, Int, Int, Int, Double, String)] -> [(Int, Int, Int, Int, Double, String)] -> Bool
yoloAcceptableDifference ref out = all (isJust . (searchForMatch out)) ref && all (isJust . (searchForMatch ref)) out
  where
    searchForMatch searchSpace (li, ri, ti, bi, p, x)
      | null matches = Nothing
      | otherwise    = Just $ head matches
      where
        -- Find all entries with the same name
        sameName = filter (\(_, _, _, _, _, x') -> x == x') searchSpace
        -- Find entries with sufficient overlap
        overlap = \(li', ri', ti', bi', _, _) -> sufficientOverlap (li, ri, ti, bi) (li', ri', ti', bi')
        overlapping = filter overlap sameName
        -- Finally, ensure that the probabilities are within a threshold
        matches = filter (\(_, _, _, _, p', _) -> abs p - p' <= 0.1) overlapping

    sufficientOverlap (li, ri, ti, bi) (li', ri', ti', bi') = overlap >= 0.6
      where
        (l, r, t, b) = (fromIntegral li, fromIntegral ri, fromIntegral ti, fromIntegral bi)
        (l', r', t', b') = (fromIntegral li', fromIntegral ri', fromIntegral ti', fromIntegral bi')
        xOverlap = maximum [ 0, minimum [r, r'] - maximum [l, l'] ]
        yOverlap = maximum [ 0, minimum [b, b'] - maximum [t, t'] ]
        overlap1 = (xOverlap * yOverlap) / ((r - l) * (b - t))
        overlap2 = (xOverlap * yOverlap) / ((r' - l') * (b' - t'))
        -- Get the least overlap, to avoid the case where a very large/small
        -- and incorrect bounding box output is treated as being correct
        overlap = minimum [overlap1, overlap2] :: Double

-- Test that the bundled TinyYOLOv2 network, stored as an ONNX file
-- loads correctly and gives expected output when given a test image
prop_yolo_loads_and_finds_objects :: Property
prop_yolo_loads_and_finds_objects = withTests 1 $ property $ do
  img <- evalIO (loadSerializedImage (imagesDir </> "person") :: IO (Maybe (S ('D3 416 416 3))))
  netPath <- evalIO $ getPathForNetwork TinyYoloV2
  net     <- evalIO $ loadTinyYoloV2 netPath

  assert $ isJust img
  assert $ isRight net

  let Just img'  = img
      Right net' = net
      probs = processOutput (runNet net' img') 0.3 2.0
      -- Values taken from running the TinyYolov2 network bundled
      -- with grenade through the ONNX inference runtime Python
      -- module
      reference = [ (75, 163, 97, 371,   0.7742048036087853, "person")
                  , (2, 86, 262, 354,    0.7187669842556927, "sheep" )
                  , (283, 409, 151, 335, 0.6768828930128936, "sheep" )
                  , (305, 410, 155, 323, 0.5562674827033209, "sheep" ) ]
  assert $ yoloAcceptableDifference reference probs

-- Test that the bundled ResNet-18 network, stored as an ONNX file,
-- loads correctly and gives expected output when given a test image
prop_resnet_loads_and_finds_objects :: Property
prop_resnet_loads_and_finds_objects = withTests 1 $ property $ do
  img <- evalIO (loadSerializedImage (imagesDir </> "dog") :: IO (Maybe (S ('D3 224 224 3))))
  netPath <- evalIO $ getPathForNetwork ResNet18
  net     <- evalIO $ loadResNet netPath

  assert $ isJust img
  assert $ isRight net

  let Just img'  = img
      Right net' = net
      S1D y = runNet net' img'
      tops  = getTop 5 $ LA.toList $ H.extract y
  tops === [162,166,167,215,164] -- beagle, Walker hound, English foxhound, Brittany spaniel, bluetick
  where
    getTop :: Ord a => Int -> [a] -> [Int]
    getTop n xs = map fst $ take n $ sortBy (flip compare `on` snd) $ zip [0..] xs

prop_superresolution_loads_and_upsamples :: Property
prop_superresolution_loads_and_upsamples = withTests 1 $ property $ do
  imgY0 <- evalIO (loadSerializedChannel (imagesDir </> "cheetah_Y0") :: IO (Maybe (S ('D3 224 224 1))))
  netPath <- evalIO $ getPathForNetwork SuperResolution
  net     <- evalIO $ loadSuperResolution netPath

  assert $ isJust imgY0
  assert $ isRight net

  let Just input = imgY0
      Right net' = net
      (S3D input') = input
      (S3D v)      = runNet net' input :: S ('D3 672 672 1)
      down         = downsample v
      diffs        = zipWith (\x y -> abs (x - y)) (concat $ D.toLists $ H.extract down) (concat $ D.toLists $ H.extract input')
      avgDiff      = (sum diffs) / (224 * 224)
  isWithinOf 0.01 avgDiff 0
  where
    downsample :: L 672 672 -> L 224 224
    downsample mat = H.build builder
      where
        mat' = H.extract mat
        builder i j = (sum vals) / 9
          where
            ro = (floor i) * 3
            co = (floor j) * 3
            vals = [mat' `D.atIndex` (ro', co') | ro' <- [ro..ro + 2], co' <- [co..co + 2]]

tests :: IO Bool
tests = checkParallel $$(discover)

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
-- | Vector is an alias to a DataFrame with order 1.
module Numeric.Vector
    (  -- * Type aliases
      Vector
    , Vec2f, Vec3f, Vec4f, Vec2d, Vec3d, Vec4d
    , Vec2i, Vec3i, Vec4i, Vec2w, Vec3w, Vec4w
      -- * Vector constructors
    , Vector2 (..), Vector3 (..), Vector4 (..)
    , DataFrame(Vec2, Vec3, Vec4)
      -- * Common operations
    , (.*.), dot, (·)
    , normL1, normL2, normLPInf, normLNInf, normLP
    , normalized
    , det2, cross, (×)
    ) where

import Numeric.DataFrame.Internal.Backend ()
import Numeric.Vector.Internal

max_stars_repo_name: achirkin/fasttensor

max_issues_repo_name: achirkin/fasttensor

max_forks_repo_name: achirkin/fasttensor