text:
{-# LANGUAGE FlexibleInstances #-}
{-|
Module      : Algorithm.ShapeConstraint
Description : evaluation of shape-constraints for symbolic regression
Copyright   : (c) Fabricio Olivetti de Franca, 2022
License     : GPL-3
Maintainer  : fabricio.olivetti@gmail.com
Stability   : experimental
Portability : POSIX

This package provides support functions to evaluate different
shape-constraints for symbolic regression using modal arithmetic.

See https://direct.mit.edu/evco/article/30/1/75/99840 for more details.

@
@article{kronberger2022shape,
  title={Shape-Constrained Symbolic Regression—Improving Extrapolation with Prior Knowledge},
  author={Kronberger, Gabriel and de Fran{\c{c}}a, Fabricio Olivetti and Burlacu, Bogdan and Haider, Christian and Kommenda, Michael},
  journal={Evolutionary computation},
  volume={30},
  number={1},
  pages={75--98},
  year={2022},
  publisher={MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info~…}
}
@

-}
module Algorithm.ShapeConstraint
    ( getViolationFun
    , Evaluator(..)
    , Shape(..)
    )
      where

import Numeric.ModalInterval
import Numeric.ModalInterval.Algorithms
import Data.SRTree

import Data.Map.Strict (Map(..))
import qualified Data.Map.Strict as M
import qualified Data.Vector as V
import Data.Bool (bool)
import Data.Maybe (fromMaybe)
import qualified Numeric.LinearAlgebra.Data as LA

-- | Data type describing the type of constraint:
--
-- * Range of the function
-- * Range of a partial derivative
-- * Monotonicity (non-increasing and non-decreasing)
-- * Partial monotonicity
-- * Inflection, Convex or Concave
--
-- The monotonicity constraints require the index of the
-- variable, the ranges require a tuple of the minimum and
-- maximum values, the partial monotonicity also requires
-- the range of the domain that has this property,
-- inflection, convexity and concavity requires the first
-- and second index of the partial derivatives.
data Shape   = Range (Double, Double)                     -- ^ f(x) \in [a, b]
             | DiffRng Int (Double, Double)               -- ^ d f(x) / dx \in [a, b]
             | NonIncreasing Int                          -- ^ d f(x) / dx \in [-inf, 0]
             | NonDecreasing Int                          -- ^ d f(x) / dx \in [0, inf]
             | PartialNonIncreasing Int (Double, Double)  -- ^ d f(x) / dx \in [-inf, 0], a <= x <= b
             | PartialNonDecreasing Int (Double, Double)  -- ^ d f(x) / dx \in [0, inf], a <= x <= b
             | Inflection Int Int                         -- ^ d^2 f(x)/dx^2 == 0
             | Convex Int Int                             -- ^ d^2 f(x)/dx^2 \in (0, Infinity)
             | Concave Int Int                            -- ^ d^2 f(x)/dx^2 \in (-Infinitiy,0)
                    deriving (Show, Read)

-- inner (a, b) no sharper than that
-- outer (a, b) no wider than that

-- | The `Evaluator` indicates the algorithm used to calculate the constraint:
--
-- * InnerInterval can result in false positives (i.e., it says it is feasible but it is not)
-- * OuterInterval can result in false negatives (i.e., it says it is infeasible but it is not)
-- * Sampling can result in false positives and usually requires a lot of samples
-- * Hybrid - not yet implemented
-- * Bisection can improve the results of OuterInterval
data Evaluator = InnerInterval | OuterInterval | Sampling Int | Hybrid | Bisection Int
                   deriving (Show, Read)

-- | Function that calculates how much of the constraints a symbolic tree vioaltes
type ConstraintFun = SRTree Int Double -> Double

-- | A list of domains of the variables
type Domains       = [(Double, Double)]

-- | A list of domains using Modal Arithmetic
type KDomains      = Map (Kaucher Double)

-- | Calculates the partial derivatives of a tree when needed.
ofShape :: (Floating a, Eq a, OptIntPow a) => SRTree Int a -> Shape -> SRTree Int a
t `ofShape` Range _ = t
t `ofShape` DiffRng ix _ = simplify $ deriveBy ix t
t `ofShape` NonIncreasing ix = simplify $ deriveBy ix t
t `ofShape` NonDecreasing ix = simplify $ deriveBy ix t
t `ofShape` PartialNonIncreasing ix _ = simplify $ deriveBy ix t
t `ofShape` PartialNonDecreasing ix _ = simplify $ deriveBy ix t
t `ofShape` Inflection ix iy = simplify $ deriveBy iy $ simplify $ deriveBy ix t
t `ofShape` Convex ix iy = simplify $ deriveBy iy $ simplify $ deriveBy ix t
t `ofShape` Concave ix iy = simplify $ deriveBy iy $ simplify $ deriveBy ix t
{-# INLINE ofShape #-}

toTuple :: Kaucher Double -> (Double, Double)
toTuple k = (fromMaybe (-1/0) $ inf k, fromMaybe (1/0) $ sup k)
{-# INLINE toTuple #-}

eps :: Double
eps = 1e-6
{-# INLINE eps #-}

calcPenalty :: Bool -> Double -> Double
calcPenalty b pnlty = bool 0 pnlty b
{-# INLINE calcPenalty #-}

inRng (a, b) (x, y) = calcPenalty (x < a - eps) (a - x) + calcPenalty (y > b + eps) (y - b)
{-# INLINE inRng #-}
isPositive (x, y)   = calcPenalty (x <= -eps) (abs x) + calcPenalty (y <= -eps) (abs y)
{-# INLINE isPositive #-}
isNegative (x, y)   = calcPenalty (x >= eps) (abs x) + calcPenalty (y >= eps) (abs y)
{-# INLINE isNegative #-}
isZero (x, y)       = calcPenalty (x <= -eps) (abs x) + calcPenalty (y >= eps) y
{-# INLINE isZero #-}

-- | Evaluates the constraint
evalConstraints :: Shape -> (Double, Double) -> Double
evalConstraints (Range ab) xy                 = inRng ab xy
evalConstraints (DiffRng _ ab) xy             = inRng ab xy
evalConstraints (NonIncreasing _) xy          = isNegative xy
evalConstraints (NonDecreasing _) xy          = isPositive xy
evalConstraints (PartialNonIncreasing _ _) xy = isNegative xy
evalConstraints (PartialNonDecreasing _ _) xy = isPositive xy
evalConstraints (Inflection _ _) xy           = isZero xy
evalConstraints (Convex _ _) xy               = isPositive xy
evalConstraints (Concave _ _) xy              = isNegative xy
{-# INLINE evalConstraints #-}

-- | Convert the domains for the partial monotonicity constraint
convertDomains :: Shape -> Map Int (Kaucher Double) -> Map Int (Kaucher Double)
convertDomains (PartialNonIncreasing ix (a, b)) ds = M.insert ix (a <.< b) ds
convertDomains (PartialNonDecreasing ix (a, b)) ds = M.insert ix (a <.< b) ds
convertDomains _ ds = ds
{-# INLINE convertDomains #-}

-- | Sum a list of constraint violations
sumCnstr :: (SRTree Int (Kaucher Double) -> Map Int (Kaucher Double) -> Kaucher Double) -> Map Int (Kaucher Double) -> [SRTree Int (Kaucher Double)] -> [Shape] -> Double -> Double
sumCnstr f ds [] []         acc = acc
sumCnstr f ds (x:xs) (z:zs) acc = sumCnstr f ds xs zs (acc + cnstr)
   where
     cnstr = evalConstraints z $ toTuple $ f x $ convertDomains z ds
{-# INLINE sumCnstr #-}

toMap :: Domains -> Map Int (Kaucher Double)
toMap    = M.fromList . zip [0..] . map (uncurry (<.<))
{-# INLINE toMap #-}

-- | Returns a constraint function given an `Evaluator` a list of `Shape` and the
-- `Domains` of the input variables.
getViolationFun :: Evaluator -> [Shape] -> Domains -> ConstraintFun
getViolationFun InnerInterval shapes domains t = sumCnstr innerApprox ds ts shapes 0.0
  where
    t'       = fmap singleton t
    ts       = map (t' `ofShape`) shapes
    ds       = toMap domains

getViolationFun OuterInterval shapes domains t = sumCnstr outerApprox ds ts shapes 0.0
  where
    t'       = fmap singleton t
    ts       = map (t' `ofShape`) shapes
    ds       = toMap domains

getViolationFun (Sampling nSamples) shapes domains t = go shapes ts samples 0.0
  where
    ds        = toMap domains
    samples   = map (M.fromList . zip [0..] . makeSamples nSamples . (`convertDomains` ds)) shapes
    getSize s = LA.size $ s M.! 0
    ts        = map (t `ofShape`) shapes

    evalSamples :: Maybe (LA.Vector Double) -> (Double, Double)
    evalSamples mxs = (fromMaybe (-1/0) mmin, fromMaybe (1/0) mmax)
      where
        mmax = LA.maxElement <$> mxs
        mmin = LA.minElement <$> mxs

    go [] _ _               acc = acc
    go (x:xs) (y:ys) (z:zs) acc = go xs ys zs (acc + cnstr)
      where
        cnstr = evalConstraints x rng
        y'    = fmap (LA.fromList . replicate (getSize z)) y
        rng   = evalSamples $ evalTreeWithMap y' z

getViolationFun e shapes domains t = error $ show e <> " evaluator not implemented"

instance OptIntPow (LA.Vector Double) where
    (^.) = (^^)
    {-# INLINE (^.) #-}

-- | Creates the samples of the `Sampling` approach
makeSamples :: Int -> Map Int (Kaucher Double) -> [LA.Vector Double]
makeSamples n domains = LA.toColumns $ LA.fromLists $ map (map midpoint) $ head $ dropWhile (\x -> length x < n) samples
  where
    step x  = let (a,b) = bisect x in [a,b]
    f       = map (concatMap step)
    samples = map sequence $ iterate f $ map (pure.snd) $ M.toAscList domains
{-# INLINE makeSamples #-}

max_stars_repo_name: folivetti/shape-constraint

max_issues_repo_name: folivetti/shape-constraint

max_forks_repo_name: folivetti/shape-constraint


text:
-- |
-- Module      : Hieroglyphs
-- Description : Cairo experiments (graphics library)
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental
-- Portability : POSIX (not sure)
--
-- Jonatan H Sundqvist
-- May 31 2015
--

-- TODO | - Lenses
--        - Performance profiling (feels chuggish sometimes, consumes a lot of CPU power)
--        - Figure out how to deal with polymorphic numerical types
--        - Polymorphic wrappers around core rendering functions
--        - This module should only deal with graphics
--        - Work out how to deal with state (logical/interaction/graphical)
--        - Wrap up graphical/animation state in single object (eg. Renderable ball with trail/colour/arrows/etc.)
--          -- Toggle visibility (eg. with checkbuttons)
--          -- Drawing bounding boxes, vector arrows
--          -- Create Drawable type
--
--        - Effects
--          -- Fading
--
--        - Worry about clean-up or rely on GC (eg. freeing surfaces) (?)
--
--        - Shape API
--          -- Should render functions simply set a path or fill them in/draw them too
--          -- Utility functions for colour/fill/stroke/line width/etc. (?)

-- SPEC | -
--        -

-- module Hieroglyphs where

---------------------------------------------------------------------------------------------------
-- We'll need these
---------------------------------------------------------------------------------------------------
import Graphics.UI.Gtk -- hiding (fill)
import qualified Graphics.Rendering.Cairo as C

import Data.Complex
import Control.Monad (when, forM_)
import Data.IORef

import qualified Copernicus as Cop
import qualified Palette

---------------------------------------------------------------------------------------------------
-- Types
---------------------------------------------------------------------------------------------------
-- TODO: One trail per body ([[Complex Double]]) (?)
data World = World { frame :: Int, size :: (Int, Int), bodies :: [Cop.Body Double], trails :: [[Complex Double]] } deriving Show

---------------------------------------------------------------------------------------------------
-- Data
---------------------------------------------------------------------------------------------------
animate = True --
π       = pi   -- :: Double -- TODO: Polymorphic (?)
τ       = 2*π  -- :: Double -- TODO: Polymorphic (?)
fps     = 30   :: Int

v = 2.0:+0.5 --40.0 :+ 20.0
g = 0.0:+(-9.82) --0.0 :+ (-9.82)

bodies' = map (\ (p', v', g') -> Cop.Body p' v' g') [(0.0:+1.0, v, g), (1.0:+0.0, (2.0):+1.50, g), (3.0:+(2.80), v, g), (2.5:+1.92, v, g)]

---------------------------------------------------------------------------------------------------
-- Functions
---------------------------------------------------------------------------------------------------
-- |
mainGTK :: IO ()
mainGTK = do
    initGUI
    window <- windowNew

    frame  <- frameNew
    canvas <- drawingAreaNew
    containerAdd frame canvas

    set window [ containerChild := frame]
    windowSetDefaultSize window 650 650

    widgetAddEvents canvas [PointerMotionMask]

    widgetShowAll window

    -- Assets
    -- planets <- C.withImageSurfaceFromPNG "assets/planets.png" $ \surface -> surface
    planets <- C.liftIO $ C.imageSurfaceCreateFromPNG "assets/planets.png"

    -- Animation
    (w,h)    <- widgetSize window
    worldVar <- newIORef $ World { frame=0, size=(w,h), bodies=bodies', trails=map (const []) bodies' } --
    when animate $ timeoutAdd (onanimate worldVar canvas) (1000 `div` fps) >> return ()

    -- Events
    canvas `on` draw $ (C.liftIO $ readIORef worldVar) >>= flip render planets -- readIORef worldVar >>= \ w -> render (fromIntegral w) (fromIntegral h) w
    canvas `on` motionNotifyEvent $ onmousemove

    window `on` configureEvent $ onresize window worldVar
    window `on` deleteEvent $ C.liftIO mainQuit >> return False -- TODO: Uhmmm... what?
    window `on` keyPressEvent $ onkeypress

    mainGUI

-- |
-- TODO: Make polymorphic
widgetSize :: WidgetClass self => self -> IO (Int, Int)
widgetSize widget = do
    w <- widgetGetAllocatedWidth widget
    h <- widgetGetAllocatedHeight widget
    return (w, h)

-- |
-- onresize :: WidgetClass self => self -> IORef World -> IO Bool
onresize window worldVar = do
    sz <- eventSize
    C.liftIO $ modifyIORef worldVar $ \ wd -> wd { size=sz }
    return True

-- |
-- onkeypress ::
onkeypress = do
    key <- eventKeyName
    C.liftIO $ print key
    return True

-- |
-- onmousemove ::
onmousemove = do
    when False $ do
        (mx,my) <- eventCoordinates
        C.liftIO $ print (mx,my)
    return False

-- Animation --------------------------------------------------------------------------------------
-- |
onanimate :: IORef World -> DrawingArea -> IO Bool
onanimate world canvas = do
    -- withMVar var print
    modifyIORef world update
    widgetQueueDraw canvas
    return True

-- |
-- TODO: Refactor, make readable (State monad, lenses)
update :: World -> World
update w@(World { frame=f, bodies=bodies'' } ) = w { frame=f+1,
                                                                       bodies=map (Cop.animate (1.0/fromIntegral fps)) bodies'',
                                                                       trails=map (\ (trail, (Cop.Body p _ _)) -> take 50 $ p:trail) $ zip (trails w) bodies'' }

-- | The time elapsed (in seconds)
-- TODO: Polymorphic (?)
elapsed :: World -> Double
elapsed world = (fromIntegral $ frame world) * 1.0/fromIntegral fps

-- Rendering --------------------------------------------------------------------------------------
-- |
-- TODO: Make polymorphic
-- TODO: Apply coordinate transformations (to World perhaps; make world a monad?)
renderBody :: Cop.Body Double -> C.Render ()
renderBody (Cop.Body (x:+y) (vx:+vy) (ax:+ay)) = do
    -- C.arc (x/3 + 200) (y/3 + 200) 12 0 τ
    -- Body
    C.arc x y 12 0 τ
    C.setSourceRGBA 0.5 0.2 0.3 1.0
    C.fill

    -- Vector arrows
    C.setSourceRGBA 0.12 0.05 1 0.8
    renderArrow (x:+y) (x:+(y+vy*10)) (abs vy*10*0.8) 5 12

    C.setSourceRGBA 1 0.05 0.12 0.8
    renderArrow (x:+y) ((x+vx*10):+y) (abs vx*10*0.8) 5 12

-- |
renderWorld :: World -> C.Render ()
renderWorld world = do

    -- X-axis
    vectorise C.moveTo $ toScreenCoords world ((-5):+0)
    vectorise C.lineTo $ toScreenCoords world (  5:+0)
    C.setSourceRGBA 1 0 0 0.8
    C.stroke

    -- Y-axis
    vectorise C.moveTo $ toScreenCoords world (0:+(-5))
    vectorise C.lineTo $ toScreenCoords world (0:+  5)
    C.setSourceRGBA 0 1 0 0.8
    C.stroke

    -- Trail(s)
    -- TODO: Better way of doing 2D 'loops'
    forM_ (zip [Palette.green, Palette.red, Palette.blue] $ trails world) $ \(fill, trail) -> do
        forM_ trail $ \ dot -> do
            -- C.setSourceRGBA (cos $ 1.3*c) (sin $ 5*c) 0.2 1.0
            Palette.choose fill
            renderCircle (toScreenCoords world dot) 3

    -- Render bodies
    forM_ (bodies world) $ renderBody . bodyToScreenCoords world

-- |
render :: World -> C.Surface -> C.Render ()
render world planets = do

    --
    when False $ do
        C.moveTo 16 44
        C.liftIO $ C.fontOptionsCreate >>= flip C.fontOptionsSetAntialias C.AntialiasSubpixel
        C.selectFontFace "Helvetica" C.FontSlantNormal C.FontWeightNormal
        C.setFontSize 30
        C.setSourceRGBA 0.62 0.62 0.62 0.7
        C.showText "Copernicus"

        extents <- C.textExtents "Copernicus"

        C.setFontSize 14
        C.setSourceRGBA 0.3 0.3 0.3 1.0
        C.moveTo 16 $ 44 + C.textExtentsHeight extents + 5
        C.showText "Jonatan H Sundqvist 2015"

    --
    renderGrid 10 10 $ fromIntegral (fst $ size world) / 10

    -- Image
    when False $ do
        {-action <- C.liftIO . C.withImageSurfaceFromPNG "assets/planets.png" $ (\surface -> return $
            do C.setSourceSurface surface 20 20 >> C.paint)-}
        -- planets <- C.liftIO $ C.imageSurfaceCreateFromPNG "assets/planets.png"
        C.setSourceSurface planets 50 50
        C.rectangle 320 310 250 250
        C.clip
        C.paint
        C.resetClip

    -- Arrow
    -- renderArrow (40:+40) (200:+120) 80 20 40
    -- TODO: Don't hard-code arrow values
    when False $ do
        forM_ (bodies world) $ \ body -> do
            let (from, to, len, width, headWidth) = (40:+50, let (Cop.Body p _ _) = body in toScreenCoords world p, 0.78 * magnitude (to-from), 40, 88)
            let thearrow = arrow from to len width headWidth
            renderArrow from to len width headWidth

            -- Arrow vertices
            C.setSourceRGBA 0 0 0 1.0
            forM_ thearrow (flip renderCircle 5)

            C.setSourceRGBA 0 0 0 1.0
            forM_ (zip [1..] thearrow) $ \(n, p) -> vectorise C.moveTo (p+(6:+6)) >> C.showText (show n)

            -- Arrow symmetry line
            vectorise C.moveTo from
            vectorise C.lineTo to
            C.setSourceRGBA 0 0 0 1.0
            C.stroke

    -- Shifting polygons
    when False $ do
        let sides       = 3 + (flip mod 10 . flip div fps $ frame world)
            radius mini = (+mini) . (*35.0) . (+1.0) . sin  $ elapsed world
            origin      = (w/2):+(h/2)
            fill        = False in forM_ [(10, (0.0, 0.5, 0, 1.0)),
                                          (30, (0.3, 0.0, 0.8, 1.0)),
                                          (50, (0.0, 0.7, 0.2, 1.0)),
                                          (70, (0.8, 0.8, 0.1, 1.0))] $ \(mini, colour) -> renderPolygon sides (radius mini) origin colour fill

    --  Bouncing balls
    renderWorld world

    -- Circle carousel
    when False $ do
        renderCircleArc 10 origin spread radius begin τ
    where count  = 10
          origin = (w/2):+(h/2)
          spread = 50 + 50 * (1 + sin (τ * rpm * elapsed world)) -- 134      -- Radius of the big circle (pixels?)
          radius = 20       -- Radius of a small circle (pixels?)
          (w,h)  = let (w', h') = size world in (fromIntegral w', fromIntegral h')
          rpm    = 0.3      --
          begin  = τ * rpm * elapsed world

-- |
renderGrid :: Int -> Int -> Double -> C.Render ()
renderGrid cols rows size = do
    sequence_ [ tilePath cl rw >> C.fill   | cl <- [0..(cols-1)], rw <- [0..(rows-1)] ] -- Tiles
    sequence_ [ tilePath cl rw >> C.stroke | cl <- [0..(cols-1)], rw <- [0..(rows-1)] ] -- Borders
    -- TODO: Figure out how to use fill AND stroke
    where chooseColour cl rw = if (cl `mod` 2) == (rw `mod` 2) then 0.3 else 0.75 -- TODO: This should be a utility function
          tilePath cl rw     = C.rectangle (fromIntegral cl*size) (fromIntegral rw*size) size size >> C.setSourceRGBA 0.22 0.81 (chooseColour cl rw) 0.32

-- |
renderArrow :: Complex Double -> Complex Double -> Double -> Double -> Double -> C.Render ()
renderArrow from to sl sw hw = do
    let (first:rest) = closePath $ arrow from to sl sw hw
    vectorise C.moveTo first
    forM_ rest $ vectorise C.lineTo

    -- C.setSourceRGBA 1 0.8 0.15 1.0
    C.setLineWidth 3
    C.stroke

-- |
-- TODO: Add arguments for colour, stroke, etc.
-- TODO: Make polymorphic
renderPolygon :: Int -> Double -> Complex Double -> (Double, Double, Double, Double) -> Bool -> C.Render ()
renderPolygon sides radius origin (r,g,b,a) filled = do
    -- TODO: Refine 'wrap-around logic'
    C.moveTo fx fy
    forM_ rest $ \(x:+y) -> C.lineTo x y

    C.setSourceRGBA r g b a
    C.setLineWidth 12
    if filled
        then C.fill
        else C.stroke
    where ((fx:+fy):rest) = polygon sides radius origin ++ [fx:+fy]

-- |
-- TODO: Options for fill/stroke, colour, width, etc.
renderCircle :: Complex Double -> Double -> C.Render ()
renderCircle (cx:+cy) radius = do
    C.arc cx cy radius 0 τ
    C.fill

-- |
-- TODO: Options for colour, width, closed/open, etc.
renderPath :: [Complex Double] -> C.Render ()
renderPath (p:ath) = do
    vectorise C.moveTo p
    forM_ ath $ vectorise C.lineTo
    C.stroke

-- |
-- Ugh, I hate underscores so much
-- TODO: Make polymorphic
renderCircleArc :: Int -> Complex Double -> Double -> Double -> Double -> Double -> C.Render ()
renderCircleArc
    count    -- Number of small circles
    (ox:+oy) -- Centre of the 'arc' (pixels?)
    spread   -- Radius of the 'arc'
    radius   -- Radius of the small circles
    begin    -- Start angle of the arc
    extent = forM_ [1..count] $ \ n -> do
        let n' = fromIntegral n
        let θ  = begin + n'*extent/fromIntegral count
        C.arc (ox - spread*cos θ) (oy - spread*sin θ) radius 0 τ
        C.setSourceRGBA (0.5 * (1 + sin θ)) (0.1*n') (1/n') 0.95
        C.fill

-- Geometry ---------------------------------------------------------------------------------------
-- |
-- TODO: Start angle
-- TODO: Invalid arguments (eg. sides < 3) (use Maybe?)
-- TODO: Make polymorphic
-- TODO: Move to geometry section
-- polygon :: (Floating f, RealFloat f) => Int -> f -> Complex f -> [Complex f]
polygon :: Int -> Double -> Complex Double -> [Complex Double]
polygon sides radius origin = [ let θ = arg n in origin + ((radius * cos θ):+(radius * sin θ)) | n <- [1..sides]]
    where arg n = fromIntegral n * (2*π)/fromIntegral sides

-- |
-- TODO: Simplify and expound comments
arrow :: Complex Double -> Complex Double -> Double -> Double -> Double -> [Complex Double]
arrow from to sl sw hw = [from     + straight (sw/2), --
                          shaftEnd + straight (sw/2), --
                          shaftEnd + straight (hw/2), --
                          to,                         --
                          shaftEnd - straight (hw/2), --
                          shaftEnd - straight (sw/2), --
                          from     - straight (sw/2)]  --
    where along a b distance = a + (mkPolar distance . snd . polar $ b-a) -- Walk distance along the from a to b
          normal a b = let (mag, arg) = polar (b-a) in mkPolar mag (arg+π/2)
          shaftEnd = along from to sl --
          straight = along (0:+0) (normal from to) -- Vector perpendicular to the centre line
          -- shaftWidth = along (0:+0) (normal a b) (sw/2) -- Half of shaft width (vector relative to symmetry line)
          -- headWidth  = along (0:+0) (normal a b) (hw/2)

-- Utilities --------------------------------------------------------------------------------------
-- |
vectorise :: (Double -> Double -> a) -> Complex Double -> a
vectorise f (re:+im) = f re im

-- |
-- TODO: Unsafe, use Maybe (?)
closePath :: [Complex Double] -> [Complex Double]
closePath path = path ++ [head path]

-- Transforms a vector from one coordinate space to another
-- by applying the given scaling and translation
-- Useful for converting between simulation and screen coordinates
-- TODO: Types for representing coordinate systems
-- TODO: Use matrix instead (scale and translate should be applied in reverse order when undoing the transformation)
-- TODO: Add undo parameter (?)
-- TODO: Pure function, move to Copernicus.hs or use Cairo matrix functions
transform :: Complex Double -> Complex Double -> Complex Double -> Complex Double
transform scaling translation = scale scaling . translate translation

-- |
scale :: Complex Double -> Complex Double -> Complex Double
scale (sx:+sy) (x:+y) = (sx*x):+(sy*y)

-- |
translate :: Complex Double -> Complex Double -> Complex Double
translate = (+)

-- |
toScreenCoords :: World -> Complex Double -> Complex Double
toScreenCoords world = transform (sx:+sy) translation
    where (w, h)   = let (w', h') = size world in (fromIntegral w', fromIntegral h')
          (sx:+sy)    = (w/10):+(-h/10)
          translation = w/(2*sx):+h/(2*sy)

-- | TODO: Implement :/
-- TODO: Dotwise application for Complex numbers
-- TODO: Don't hard-code the world-size of the canvas (currently 10-by-10)
toWorldCoords :: World -> Complex Double -> Complex Double
toWorldCoords world = translate (-translation) . scale (sx**(-1):+sy**(-1))
    where (w, h) = let (w', h') = size world in (fromIntegral w', fromIntegral h')
          (sx:+sy)    = (w/10):+(-h/10)
          translation = (w/(2*sx)):+(h/(2*sy))

-- |
bodyToScreenCoords :: World -> Cop.Body Double -> Cop.Body Double
bodyToScreenCoords world (Cop.Body p v' a') = Cop.Body (toScreenCoords world p) v' a' -- TODO: Make transformation less ad-hoc

---------------------------------------------------------------------------------------------------
-- Entry point
---------------------------------------------------------------------------------------------------
main :: IO ()
main = mainGTK

ext: hs

max_stars_repo_name: SwiftsNamesake/Copernicus

max_issues_repo_name: SwiftsNamesake/Copernicus

max_forks_repo_name: SwiftsNamesake/Copernicus


text:
module FractalStream.Models
    ( Coordinate
    , ComplexParametric1d
    ) where

import Lang.Numbers

type family Coordinate model :: *
data ComplexParametric1d

{-
import Data.Word
import Data.Complex

newtype ViewCoordinate a = ViewCoordinate a

data ParameterType = N | Z | R | C

data Parameter = Parameter
  { parameterType :: !ParameterType
  , parameterName :: !String
  , parameterDesc :: !String
  }

data Group = Group
  { groupTitle  :: !String
  , groupParams :: [Parameter]
  }

epsilon :: Parameter
epsilon = Parameter
  { parameterType = R
  , parameterName = "Epsilon"
  , parameterDesc = concat
      [ "How close together two complex numbers should"
      , " be in order for them to be considered equal?" ]
  }

infinity :: Parameter
infinity = Parameter
  { parameterType = R
  , parameterName = "Infinity"
  , parameterDesc = concat
      [ "How large a number should be in order to"
      , " consider it escaped / near infinity?" ]
  }

maxIter :: Parameter
maxIter = Parameter
  { parameterType = N
  , parameterName = "Maximum iteration count"
  , parameterDesc = concat
      [ "How many iterations should we try before giving up?" ]
  }

data Int32

data ParameterDesc where
  = Param Parameter
  |

complexParametric1d
  =  viewCoordinate @(Complex Double) "C"
  <> parameter @Int32 "Maximum iteration count"
  <> parameter @Double "Epsilon"
  <> parameter @Double "Infinity"

complexParametricDynamics1d
  =  viewCoordinate @(Complex Double) "Z"
  <> parameter @(Complex Double) "C"
  <> parameter @Int32 "Maximum iteration count"
  <> parameter @Double "Epsilon"
  <> parameter @Double "Infinity"

complexDynamics1d
  =  viewCoordinate @(Complex Double) "Z"
  <> parameter @Int32 "Maximum iteration count"
  <> parameter @Double "Epsilon"
  <> parameter @Double "Infinity"
-}

type instance Coordinate ComplexParametric1d = C --(Double, Double)

max_stars_repo_name: matt-noonan/FractalStream

max_issues_repo_name: matt-noonan/FractalStream

max_forks_repo_name: matt-noonan/FractalStream


text:
{-# LANGUAGE CPP                                      #-}
{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE PolyKinds                                #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE ViewPatterns                             #-}
{-# OPTIONS_GHC -fno-warn-orphans                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

#if MIN_VERSION_base(4,12,0)
{-# LANGUAGE NoStarIsType #-}
#endif

-- |
-- Module      : Numeric.LinearAlgebra.Static.Backprop
-- Copyright   : (c) Justin Le 2018
-- License     : BSD3
--
-- Maintainer  : justin@jle.im
-- Stability   : experimental
-- Portability : non-portable
--
-- A wrapper over "Numeric.LinearAlgebra.Static" (type-safe vector and
-- matrix operations based on blas/lapack) that allows its operations to
-- work with <https://hackage.haskell.org/package/backprop backprop>.  Also
-- provides orphan instances of 'Backprop' for types in
-- "Numeric.LinearAlgebra.Static".
--
-- In short, these functions are "lifted" to work with 'BVar's.
--
-- Using 'evalBP' will run the original operation:
--
-- @
-- 'evalBP' :: (forall s. 'Reifies' s 'W'. 'BVar' s a -> 'BVar' s b) -> a -> b
-- @
--
-- But using 'gradBP' or 'backprop' will give you the gradient:
--
-- @
-- 'gradBP' :: (forall s. 'Reifies' s 'W'. 'BVar' s a -> 'BVar' s b) -> a -> a
-- @
--
-- These can act as a drop-in replacement to the API of
-- "Numeric.LinearAlgebra.Static".  Just change your imports, and your
-- functions are automatically backpropagatable.  Useful types are all
-- re-exported.
--
-- Also contains 'sumElements' 'BVar' operation.
--
-- Formulas for gradients come from the following papers:
--
--     * https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf
--     * http://www.dtic.mil/dtic/tr/fulltext/u2/624426.pdf
--     * http://www.cs.cmu.edu/~zkolter/course/15-884/linalg-review.pdf
--     * https://arxiv.org/abs/1602.07527
--
-- Some functions are notably unlifted:
--
--     * 'H.svd': I can't find any resources that allow you to backpropagate
--     if the U and V matrices are used!  If you find one, let me know, or
--     feel free to submit a PR!  Because of this, Currently only a version
--     that exports only the singular values is exported.
--     * 'H.svdTall', 'H.svdFlat': Not sure where to start for these
--     * 'qr': Same story.
--     https://github.com/tensorflow/tensorflow/issues/6504 might yield
--     a clue?
--     * 'H.her': No 'Num' instance for 'H.Her' makes this impossible at
--     the moment with the current backprop API
--     * 'H.exmp': Definitely possible, but I haven't dug deep enough to
--     figure it out yet!  There is a description here
--     https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf but it
--     requires some things I am not familiar with yet.  Feel free to
--     submit a PR!
--     * 'H.sqrtm': Also likely possible.  Maybe try to translate
--     http://people.cs.umass.edu/~smaji/projects/matrix-sqrt/ ?  PRs
--     welcomed!
--     * 'H.linSolve': Haven't figured out where to start!
--     * 'H.</>': Same story
--     * Functions returning existential types, like 'H.withNullSpace',
--     'H.withOrth', 'H.withRows', etc.; not quite sure what the best way
--     to handle these are at the moment.
--     * 'H.withRows' and 'H.withColumns' made "type-safe", without
--     existential types, with 'fromRows' and 'fromColumns'.

module Numeric.LinearAlgebra.Static.Backprop (
  -- * Vector
    H.R
  , H.ℝ
  , vec2
  , vec3
  , vec4
  , (&)
  , (#)
  , split
  , headTail
  , vector
  , linspace
  , H.range
  , H.dim
  -- * Matrix
  , H.L
  , H.Sq
  , row
  , col
  , (|||)
  , (===)
  , splitRows
  , splitCols
  , unrow
  , uncol
  , tr
  , H.eye
  , diag
  , matrix
  -- * Complex
  , H.ℂ
  , H.C
  , H.M
  , H.𝑖
  , toComplex
  , fromComplex
  , complex
  , real
  , imag
  , sqMagnitude
  , magnitude
  -- * Products
  , (<>)
  , (#>)
  , (<.>)
  -- * Factorizations
  , svd
  , svd_
  , H.Eigen
  , eigensystem
  , eigenvalues
  , chol
  -- * Norms
  , H.Normed
  , norm_0
  , norm_1V
  , norm_1M
  , norm_2V
  , norm_2M
  , norm_InfV
  , norm_InfM
  -- * Misc
  , mean
  , meanCov
  , meanL
  , cov
  , H.Disp(..)
  -- ** Domain
  , H.Domain
  , mul
  , app
  , dot
  , cross
  , diagR
  , vmap
  , vmap'
  , dvmap
  , mmap
  , mmap'
  , dmmap
  , outer
  , zipWithVector
  , zipWithVector'
  , dzipWithVector
  , det
  , invlndet
  , lndet
  , inv
  -- ** Conversions
  , toRows
  , toColumns
  , fromRows
  , fromColumns
  -- ** Misc Operations
  , konst
  , sumElements
  , extractV
  , extractM
  , create
  , H.Diag
  , takeDiag
  , H.Sym
  , sym
  , mTm
  , unSym
  , (<·>)
  -- * Backprop types re-exported
  -- | Re-exported for convenience.
  --
  -- @since 0.1.1.0
  , BVar
  , Backprop
  , Reifies
  , W
  ) where

import           Data.Bifunctor
import           Data.Coerce
import           Data.Functor.Identity
import           Data.Maybe
import           Data.Proxy
import           Data.Vinyl                          (Rec(..))
import           Foreign.Storable
import           GHC.TypeLits
import           Lens.Micro hiding                   ((&))
import           Numeric.Backprop
import           Numeric.Backprop.Class
import           Unsafe.Coerce
import qualified Data.Vector                         as V
import qualified Data.Vector.Generic                 as VG
import qualified Data.Vector.Generic.Sized           as SVG
import qualified Data.Vector.Sized                   as SV
import qualified Data.Vector.Storable.Sized          as SVS
import qualified Numeric.Backprop.Explicit           as BE
import qualified Numeric.LinearAlgebra               as HU
import qualified Numeric.LinearAlgebra.Static        as H
import qualified Numeric.LinearAlgebra.Static.Vector as H
import qualified Prelude.Backprop                    as B

#if MIN_VERSION_base(4,11,0)
import           Prelude hiding               ((<>))
#endif

instance Backprop (H.R n) where
    zero = zeroNum
    add = addNum
    one = oneNum

instance Backprop (H.C n) where
    zero = zeroNum
    add = addNum
    one = oneNum

instance (KnownNat n, KnownNat m) => Backprop (H.L n m) where
    zero = zeroNum
    add = addNum
    one = oneNum

instance (KnownNat n, KnownNat m) => Backprop (H.M n m) where
    zero = zeroNum
    add = addNum
    one = oneNum

instance KnownNat n => Backprop (H.Sym n) where
    zero = zeroNum
    add = addNum
    one = oneNum

vec2
    :: Reifies s W
    => BVar s H.ℝ
    -> BVar s H.ℝ
    -> BVar s (H.R 2)
vec2 = isoVar2 H.vec2 (\(H.rVec->v) -> (SVS.index v 0, SVS.index v 1))
{-# INLINE vec2 #-}

vec3
    :: Reifies s W
    => BVar s H.ℝ
    -> BVar s H.ℝ
    -> BVar s H.ℝ
    -> BVar s (H.R 3)
vec3 = isoVar3 H.vec3 (\(H.rVec->v) -> (SVS.index v 0, SVS.index v 1, SVS.index v 2))
{-# INLINE vec3 #-}

vec4
    :: Reifies s W
    => BVar s H.ℝ
    -> BVar s H.ℝ
    -> BVar s H.ℝ
    -> BVar s H.ℝ
    -> BVar s (H.R 4)
vec4 vX vY vZ vW = isoVarN
    (\(Identity x :& Identity y :& Identity z :& Identity w :& RNil) -> H.vec4 x y z w)
    (\(H.rVec->v) -> Identity (SVS.index v 0)
                  :& Identity (SVS.index v 1)
                  :& Identity (SVS.index v 2)
                  :& Identity (SVS.index v 3)
                  :& RNil
    )
    (vX :& vY :& vZ :& vW :& RNil)
{-# INLINE vec4 #-}

(&) :: (KnownNat n, 1 <= n, KnownNat (n + 1), Reifies s W)
    => BVar s (H.R n)
    -> BVar s H.ℝ
    -> BVar s (H.R (n + 1))
(&) = isoVar2 (H.&) (\(H.split->(dxs,dy)) -> (dxs, fst (H.headTail dy)))
infixl 4 &
{-# INLINE (&) #-}

(#) :: (KnownNat n, KnownNat m, Reifies s W)
    => BVar s (H.R n)
    -> BVar s (H.R m)
    -> BVar s (H.R (n + m))
(#) = isoVar2 (H.#) H.split
infixl 4 #
{-# INLINE (#) #-}

split
    :: forall p n s. (KnownNat p, KnownNat n, p <= n, Reifies s W)
    => BVar s (H.R n)
    -> (BVar s (H.R p), BVar s (H.R (n - p)))
split v = (t ^^. _1, t ^^. _2)
  where
    t = isoVar H.split (uncurry (H.#)) v
    {-# NOINLINE t #-}
{-# INLINE split #-}

headTail
    :: (Reifies s W, KnownNat n, 1 <= n)
    => BVar s (H.R n)
    -> (BVar s H.ℝ, BVar s (H.R (n - 1)))
headTail v = (t ^^. _1, t ^^. _2)
  where
    t = isoVar H.headTail
                   (\(d, dx) -> (H.konst d :: H.R 1) H.# dx)
                   v
    {-# NOINLINE t #-}
{-# INLINE headTail #-}

-- | Potentially extremely bad for anything but short lists!!!
vector
    :: forall n s. (KnownNat n, Reifies s W)
    => SV.Vector n (BVar s H.ℝ)
    -> BVar s (H.R n)
vector = BE.isoVar afSV
            (H.vecR . SVG.convert) (SVG.convert . H.rVec)
       . collectVar
{-# INLINE vector #-}

linspace
    :: forall n s. (KnownNat n, Reifies s W)
    => BVar s H.ℝ
    -> BVar s H.ℝ
    -> BVar s (H.R n)
linspace = liftOp2 . op2 $ \l u ->
    ( H.linspace (l, u)
    , \d -> let n1 = fromInteger $ natVal (Proxy @n) - 1
                dDot = ((H.range - 1) H.<.> d) / n1
                dSum = HU.sumElements . H.extract $ d
            in  (dSum - dDot, dDot)
    )
{-# INLINE linspace #-}

row :: Reifies s W
    => BVar s (H.R n)
    -> BVar s (H.L 1 n)
row = isoVar H.row H.unrow
{-# INLINE row #-}

col :: (KnownNat n, Reifies s W)
    => BVar s (H.R n)
    -> BVar s (H.L n 1)
col = isoVar H.col H.uncol
{-# INLINE col #-}

(|||) :: (KnownNat c, KnownNat r1, KnownNat (r1 + r2), Reifies s W)
      => BVar s (H.L c r1)
      -> BVar s (H.L c r2)
      -> BVar s (H.L c (r1 + r2))
(|||) = isoVar2 (H.|||) H.splitCols
infixl 3 |||
{-# INLINE (|||) #-}

(===) :: (KnownNat c, KnownNat r1, KnownNat (r1 + r2), Reifies s W)
      => BVar s (H.L r1        c)
      -> BVar s (H.L r2        c)
      -> BVar s (H.L (r1 + r2) c)
(===) = isoVar2 (H.===) H.splitRows
infixl 2 ===
{-# INLINE (===) #-}

splitRows
    :: forall p m n s. (KnownNat p, KnownNat m, KnownNat n, p <= m, Reifies s W)
    => BVar s (H.L m n)
    -> (BVar s (H.L p n), BVar s (H.L (m - p) n))
splitRows v = (t ^^. _1, t ^^. _2)
  where
    t = isoVar H.splitRows (uncurry (H.===)) v
    {-# NOINLINE t #-}
{-# INLINE splitRows #-}

splitCols
    :: forall p m n s. (KnownNat p, KnownNat m, KnownNat n, KnownNat (n - p), p <= n, Reifies s W)
    => BVar s (H.L m n)
    -> (BVar s (H.L m p), BVar s (H.L m (n - p)))
splitCols v = (t ^^. _1, t ^^. _2)
  where
    t = isoVar H.splitCols (uncurry (H.|||)) v
    {-# NOINLINE t #-}
{-# INLINE splitCols #-}

unrow
    :: (KnownNat n, Reifies s W)
    => BVar s (H.L 1 n)
    -> BVar s (H.R n)
unrow = isoVar H.unrow H.row
{-# INLINE unrow #-}

uncol
    :: (KnownNat n, Reifies s W)
    => BVar s (H.L n 1)
    -> BVar s (H.R n)
uncol = isoVar H.uncol H.col
{-# INLINE uncol #-}

tr  :: (HU.Transposable m mt, HU.Transposable mt m, Backprop m, Reifies s W)
    => BVar s m
    -> BVar s mt
tr = isoVar H.tr H.tr
{-# INLINE tr #-}

diag
    :: (KnownNat n, Reifies s W)
    => BVar s (H.R n)
    -> BVar s (H.Sq n)
diag = liftOp1 . op1 $ \x -> (H.diag x, H.takeDiag)
{-# INLINE diag #-}

-- | Potentially extremely bad for anything but short lists!!!
matrix
    :: forall m n s. (KnownNat m, KnownNat n, Reifies s W)
    => [BVar s H.ℝ]
    -> BVar s (H.L m n)
matrix = maybe (error "matrix: invalid number of elements")
               ( isoVar (H.vecL . SVG.convert . runABP) (ABP . SVG.convert . H.lVec)
               . collectVar
               . ABP
               )
       . SV.fromList @(m * n)
{-# INLINE matrix #-}

-- | Matrix product
(<>)
    :: (KnownNat m, KnownNat k, KnownNat n, Reifies s W)
    => BVar s (H.L m k)
    -> BVar s (H.L k n)
    -> BVar s (H.L m n)
(<>) = mul
infixr 8 <>
{-# INLINE (<>) #-}

-- | Matrix-vector product
(#>)
    :: (KnownNat m, KnownNat n, Reifies s W)
    => BVar s (H.L m n)
    -> BVar s (H.R n)
    -> BVar s (H.R m)
(#>) = app
infixr 8 #>
{-# INLINE (#>) #-}

-- | Dot product
(<.>)
    :: (KnownNat n, Reifies s W)
    => BVar s (H.R n)
    -> BVar s (H.R n)
    -> BVar s H.ℝ
(<.>) = dot
infixr 8 <.>
{-# INLINE (<.>) #-}

-- | Can only get the singular values, for now.  Let me know if you find an
-- algorithm that can compute the gradients based on differentials for the
-- other matricies!
--
svd :: forall m n s. (KnownNat m, KnownNat n, Reifies s W)
    => BVar s (H.L m n)
    -> BVar s (H.R n)
svd = liftOp1 . op1 $ \x ->
    let (u, σ, v) = H.svd x
    in  ( σ
        , \(dΣ :: H.R n) -> (u H.<> H.diagR 0 dΣ) H.<> H.tr v
                -- must manually associate because of bug in diagR in
                -- hmatrix-0.18.2.0
        )
{-# INLINE svd #-}

-- | Version of 'svd' that returns the full SVD, but if you attempt to find
-- the gradient, it will fail at runtime if you ever use U or V.
svd_
    :: forall m n s. (KnownNat m, KnownNat n, Reifies s W)
    => BVar s (H.L m n)
    -> (BVar s (H.L m m), BVar s (H.R n), BVar s (H.L n n))
svd_ r = (t ^^. _1, t ^^. _2, t ^^. _3)
  where
    o :: Op '[H.L m n] (H.L m m, H.R n, H.L n n)
    o = op1 $ \x ->
        let msv@(u, _, v) = H.svd x
        in  ( msv
            , \(dU, dΣ, dV) ->
                    if H.norm_0 dU == 0 && H.norm_0 dV == 0
                      then (u H.<> H.diagR 0 dΣ) H.<> H.tr v
                      else error "svd_: Cannot backprop if U and V are used."
            )
    {-# INLINE o #-}
    t = liftOp1 o r
    {-# NOINLINE t #-}
{-# INLINE svd_ #-}

helpEigen :: KnownNat n => H.Sym n -> (H.R n, H.L n n, H.L n n, H.L n n)
helpEigen x = (l, v, H.inv v, H.tr v)
  where
    (l, v) = H.eigensystem x
{-# INLINE helpEigen #-}

-- | /NOTE/ The gradient is not necessarily symmetric!  The gradient is not
-- meant to be retireved directly; insteadl, 'eigenvalues' is meant to be
-- used as a part of a larger computation, and the gradient as an
-- intermediate step.
eigensystem
    :: forall n s. (KnownNat n, Reifies s W)
    => BVar s (H.Sym n)
    -> (BVar s (H.R n), BVar s (H.L n n))
eigensystem u = (t ^^. _1, t ^^. _2)
  where
    o :: Op '[H.Sym n] (H.R n, H.L n n)
    o = op1 $ \x ->
        let (l, v, vInv, vTr) = helpEigen x
            lRep = H.rowsL . SV.replicate $ l
            fMat = (1 - H.eye) * (lRep - H.tr lRep)
        in  ( (l, v)
            , \(dL, dV) -> unsafeCoerce $
                       H.tr vInv
                  H.<> (H.diag dL + fMat * (vTr H.<> dV))
                  H.<> vTr
            )
    {-# INLINE o #-}
    t = liftOp1 o u
    {-# NOINLINE t #-}
{-# INLINE eigensystem #-}

-- | /NOTE/ The gradient is not necessarily symmetric!  The gradient is not
-- meant to be retireved directly; insteadl, 'eigenvalues' is meant to be
-- used as a part of a larger computation, and the gradient as an
-- intermediate step.
eigenvalues
    :: forall n s. (KnownNat n, Reifies s W)
    => BVar s (H.Sym n)
    -> BVar s (H.R n)
eigenvalues = liftOp1 . op1 $ \x ->
    let (l, _, vInv, vTr) = helpEigen x
    in  ( l
        , \dL -> unsafeCoerce $
                 H.tr vInv H.<> H.diag dL H.<> vTr
        )
{-# INLINE eigenvalues #-}

-- | Algorithm from https://arxiv.org/abs/1602.07527
--
-- The paper also suggests a potential imperative algorithm that might
-- help.  Need to benchmark to see what is best.
--
-- /NOTE/ The gradient is not necessarily symmetric!  The gradient is not
-- meant to be retireved directly; insteadl, 'eigenvalues' is meant to be
-- used as a part of a larger computation, and the gradient as an
-- intermediate step.
chol
    :: forall n s. (KnownNat n, Reifies s W)
    => BVar s (H.Sym n)
    -> BVar s (H.Sq n)
chol = liftOp1 . op1 $ \x ->
    let l = H.chol x
        lInv = H.inv l
        phi :: H.Sq n
        phi = H.build $ \i j -> case compare i j of
                                  LT -> 1
                                  EQ -> 0.5
                                  GT -> 0
    in  ( l
        , \dL -> let s = H.tr lInv H.<> (phi * (H.tr l H.<> dL)) H.<> lInv
                 in  unsafeCoerce $ s + H.tr s - H.eye * s
        )
{-# INLINE chol #-}

-- | Number of non-zero items
norm_0
    :: (H.Normed a, Backprop a, Reifies s W)
    => BVar s a
    -> BVar s H.ℝ
norm_0 = liftOp1 . op1 $ \x -> (H.norm_0 x, const (zero x))
{-# INLINE norm_0 #-}

-- | Sum of absolute values
norm_1V
    :: (KnownNat n, Reifies s W)
    => BVar s (H.R n)
    -> BVar s H.ℝ
norm_1V = liftOp1 . op1 $ \x -> (H.norm_1 x, (* signum x) . H.konst)
{-# INLINE norm_1V #-}

-- | Maximum 'H.norm_1' of columns
norm_1M
    :: (KnownNat n, KnownNat m, Reifies s W)
    => BVar s (H.L n m)
    -> BVar s H.ℝ
norm_1M = liftOp1 . op1 $ \x ->
    let n = H.norm_1 x
    in  (n, \d -> let d' = H.konst d
                  in  H.colsL
                    . SV.map (\c -> if H.norm_1 c == n
                                      then d' * signum c
                                      else 0
                             )
                    . H.lCols
                    $ x
        )
{-# INLINE norm_1M #-}

-- | Square root of sum of squares
--
-- Be aware that gradient diverges when the norm is zero
norm_2V
    :: (KnownNat n, Reifies s W)
    => BVar s (H.R n)
    -> BVar s H.ℝ
norm_2V = liftOp1 . op1 $ \x ->
    let n = H.norm_2 x
    in (n, \d -> x * H.konst (d / n))
{-# INLINE norm_2V #-}

-- | Maximum singular value
norm_2M
    :: (KnownNat n, KnownNat m, Reifies s W)
    => BVar s (H.L n m)
    -> BVar s H.ℝ
norm_2M = liftOp1 . op1 $ \x ->
    let n = H.norm_2 x
        (head.H.toColumns->u1,_,head.H.toColumns->v1) = H.svd x
    in (n, \d -> H.konst d * (u1 `H.outer` v1))
{-# INLINE norm_2M #-}

-- | Maximum absolute value
norm_InfV
    :: (KnownNat n, Reifies s W)
    => BVar s (H.R n)
    -> BVar s H.ℝ
norm_InfV = liftOp1 . op1 $ \x ->
    let n :: H.ℝ
        n = H.norm_Inf x
    in  (n, \d -> H.vecR
                . SVS.map (\e -> if abs e == n
                                   then signum e * d
                                   else 0
                          )
                . H.rVec
                $ x
        )
{-# ANN norm_InfV "HLint: ignore Use camelCase" #-}
{-# INLINE norm_InfV #-}

-- | Maximum 'H.norm_1' of rows
norm_InfM
    :: (KnownNat n, KnownNat m, Reifies s W)
    => BVar s (H.L n m)
    -> BVar s H.ℝ
norm_InfM = liftOp1 . op1 $ \x ->
    let n = H.norm_Inf x
    in  (n, \d -> let d' = H.konst d
                  in  H.rowsL
                    . SV.map (\c -> if H.norm_1 c == n
                                      then d' * signum c
                                      else 0
                             )
                    . H.lRows
                    $ x
        )
{-# ANN norm_InfM "HLint: ignore Use camelCase" #-}
{-# INLINE norm_InfM #-}

mean
    :: (KnownNat n, 1 <= n, Reifies s W)
    => BVar s (H.R n)
    -> BVar s H.ℝ
mean = liftOp1 . op1 $ \x -> (H.mean x, H.konst . (/ H.norm_0 x))
{-# INLINE mean #-}

gradCov
    :: forall m n. (KnownNat m, KnownNat n)
    => H.L m n
    -> H.R n
    -> H.Sym n
    -> H.L m n
gradCov x μ dσ = H.rowsL
               . SV.map (subtract (dDiffsSum / m))
               . H.lRows
               $ dDiffs
  where
    diffs = H.rowsL . SV.map (subtract μ) . H.lRows $ x
    dDiffs = H.konst (2/n) * (diffs H.<> H.tr (H.unSym dσ))
    dDiffsSum = sum . H.toRows $ dDiffs
    m = fromIntegral $ natVal (Proxy @m)
    n = fromIntegral $ natVal (Proxy @n)
{-# INLINE gradCov #-}

-- | Mean and covariance.  If you know you only want to use one or the
-- other, use 'meanL' or 'cov'.
meanCov
    :: forall m n s. (KnownNat n, KnownNat m, 1 <= m, Reifies s W)
    => BVar s (H.L m n)
    -> (BVar s (H.R n), BVar s (H.Sym n))
meanCov v = (t ^^. _1, t ^^. _2)
  where
    m = fromInteger $ natVal (Proxy @m)
    t = ($ v) . liftOp1 . op1 $ \x ->
        let ms@(μ, _) = H.meanCov x
        in  ( ms
            , \(dμ, dσ) ->
                let gradMean = H.rowsL
                             . SV.replicate
                             $ (dμ / H.konst m)
                in  gradMean + gradCov x μ dσ
            )
    {-# NOINLINE t #-}
{-# INLINE meanCov #-}

-- | 'meanCov', but if you know you won't use the covariance.
meanL
    :: forall m n s. (KnownNat n, KnownNat m, 1 <= m, Reifies s W)
    => BVar s (H.L m n)
    -> BVar s (H.R n)
meanL = liftOp1 . op1 $ \x ->
    ( fst (H.meanCov x)
    , H.rowsL . SV.replicate . (/ H.konst m)
    )
  where
    m = fromInteger $ natVal (Proxy @m)
{-# INLINE meanL #-}

-- | 'cov', but if you know you won't use the covariance.
cov
    :: forall m n s. (KnownNat n, KnownNat m, 1 <= m, Reifies s W)
    => BVar s (H.L m n)
    -> BVar s (H.Sym n)
cov = liftOp1 . op1 $ \x ->
    let (μ, σ) = H.meanCov x
    in  (σ, gradCov x μ)
{-# INLINE cov #-}

mul :: ( KnownNat m
       , KnownNat k
       , KnownNat n
       , H.Domain field vec mat
       , Backprop (mat m k)
       , Backprop (mat k n)
       , HU.Transposable (mat m k) (mat k m)
       , HU.Transposable (mat k n) (mat n k)
       , Reifies s W
       )
    => BVar s (mat m k)
    -> BVar s (mat k n)
    -> BVar s (mat m n)
mul = liftOp2 . op2 $ \x y ->
    ( x `H.mul` y
    , \d -> (d `H.mul` H.tr y, H.tr x `H.mul` d)
    )
{-# INLINE mul #-}

app :: ( KnownNat m
       , KnownNat n
       , H.Domain field vec mat
       , HU.Transposable (mat m n) (mat n m)
       , Backprop (mat m n)
       , Backprop (vec n)
       , Reifies s W
       )
    => BVar s (mat m n)
    -> BVar s (vec n)
    -> BVar s (vec m)
app = liftOp2 . op2 $ \xs y ->
    ( xs `H.app` y
    , \d -> (d `H.outer` y, H.tr xs `H.app` d)
    )
{-# INLINE app #-}

dot :: ( KnownNat n
       , H.Domain field vec mat
       , H.Sized field (vec n) d
       , Num (vec n)
       , Backprop (vec n)
       , Reifies s W
       )
    => BVar s (vec n)
    -> BVar s (vec n)
    -> BVar s field
dot = liftOp2 . op2 $ \x y ->
    ( x `H.dot` y
    , \d -> let d' = H.konst d
            in  (d' * y, x * d')
    )
{-# INLINE dot #-}

cross
    :: ( H.Domain field vec mat
       , Reifies s W
       , Backprop (vec 3)
       )
    => BVar s (vec 3)
    -> BVar s (vec 3)
    -> BVar s (vec 3)
cross = liftOp2 . op2 $ \x y ->
    ( x `H.cross` y
    , \d -> (y `H.cross` d, d `H.cross` x)
    )
{-# INLINE cross #-}

-- | Create matrix with diagonal, and fill with default entries
diagR
    :: forall m n k field vec mat s.
       ( H.Domain field vec mat
       , Num (vec k)
       , Num (mat m n)
       , KnownNat m
       , KnownNat n
       , KnownNat k
       , HU.Container HU.Vector field
       , H.Sized field (mat m n) HU.Matrix
       , H.Sized field (vec k) HU.Vector
       , Backprop field
       , Backprop (vec k)
       , Reifies s W
       )
    => BVar s field             -- ^ default value
    -> BVar s (vec k)           -- ^ diagonal
    -> BVar s (mat m n)
diagR = liftOp2 . op2 $ \c x ->
    ( H.diagR c x
    , \d -> ( HU.sumElements . H.extract $ H.diagR 1 (0 :: vec k) * d
            , fromJust . H.create . HU.takeDiag . H.extract $ d
            )
    )
{-# INLINE diagR #-}

-- | Note: if possible, use the potentially much more performant 'vmap''.
vmap
    :: (KnownNat n, Reifies s W)
    => (BVar s H.ℝ -> BVar s H.ℝ)
    -> BVar s (H.R n)
    -> BVar s (H.R n)
vmap f = isoVar (H.vecR . SVG.convert @V.Vector . runABP)
                (ABP . SVG.convert . H.rVec)
       . B.fmap f
       . isoVar (ABP . SVG.convert . H.rVec) (H.vecR . SVG.convert . runABP)
{-# INLINE vmap #-}

-- | 'vmap', but potentially more performant.  Only usable if the mapped
-- function does not depend on any external 'BVar's.
vmap'
    :: ( Num (vec n)
       , Storable field
       , H.Sized field (vec n) HU.Vector
       , Backprop (vec n)
       , Backprop field
       , Reifies s W
       )
    => (forall s'. Reifies s' W => BVar s' field -> BVar s' field)
    -> BVar s (vec n)
    -> BVar s (vec n)
vmap' f = liftOp1 . op1 $ bimap (fromJust . H.create . VG.convert)
                                ((*) . fromJust . H.create . VG.convert)
                        . V.unzip
                        . V.map (backprop f)
                        . VG.convert
                        . H.extract
{-# INLINE vmap' #-}

-- TODO: Can be made more efficient if backprop exports
-- a custom-total-derivative version

-- | Note: Potentially less performant than 'vmap''.
dvmap
    :: ( KnownNat n
       , H.Domain field vec mat
       , Num (vec n)
       , Backprop (vec n)
       , Backprop field
       , Reifies s W
       )
    => (forall s'. Reifies s' W => BVar s' field -> BVar s' field)
    -> BVar s (vec n)
    -> BVar s (vec n)
dvmap f = liftOp1 . op1 $ \x ->
    ( H.dvmap (evalBP f) x
    , (H.dvmap (gradBP f) x *)
    )
{-# INLINE dvmap #-}

-- | Note: if possible, use the potentially much more performant 'mmap''.
mmap
    :: (KnownNat n, KnownNat m, Reifies s W)
    => (BVar s H.ℝ -> BVar s H.ℝ)
    -> BVar s (H.L n m)
    -> BVar s (H.L n m)
mmap f = isoVar (H.vecL . SVG.convert @V.Vector . runABP)
                (ABP . SVG.convert . H.lVec)
       . B.fmap f
       . isoVar (ABP . SVG.convert . H.lVec) (H.vecL . SVG.convert . runABP)
{-# INLINE mmap #-}

-- | 'mmap', but potentially more performant.  Only usable if the mapped
-- function does not depend on any external 'BVar's.
mmap'
    :: forall n m mat field s.
       ( KnownNat m
       , Num (mat n m)
       , Backprop (mat n m)
       , Backprop field
       , H.Sized field (mat n m) HU.Matrix
       , HU.Element field
       , Reifies s W
       )
    => (forall s'. Reifies s' W => BVar s' field -> BVar s' field)
    -> BVar s (mat n m)
    -> BVar s (mat n m)
mmap' f = liftOp1 . op1 $ bimap (fromJust . H.create . HU.reshape m . VG.convert)
                                ((*) . fromJust . H.create . HU.reshape m . VG.convert)
                        . V.unzip
                        . V.map (backprop f)
                        . VG.convert
                        . HU.flatten
                        . H.extract
  where
    m :: Int
    m = fromInteger $ natVal (Proxy @m)
{-# INLINE mmap' #-}

-- | Note: Potentially less performant than 'mmap''.
dmmap
    :: ( KnownNat n
       , KnownNat m
       , H.Domain field vec mat
       , Num (mat n m)
       , Backprop (mat n m)
       , Backprop field
       , Reifies s W
       )
    => (forall s'. Reifies s' W => BVar s' field -> BVar s' field)
    -> BVar s (mat n m)
    -> BVar s (mat n m)
dmmap f = liftOp1 . op1 $ \x ->
    ( H.dmmap (evalBP f) x
    , (H.dmmap (gradBP f) x *)
    )
{-# INLINE dmmap #-}

outer
    :: ( KnownNat m
       , KnownNat n
       , H.Domain field vec mat
       , HU.Transposable (mat n m) (mat m n)
       , Backprop (vec n)
       , Backprop (vec m)
       , Reifies s W
       )
    => BVar s (vec n)
    -> BVar s (vec m)
    -> BVar s (mat n m)
outer = liftOp2 . op2 $ \x y ->
    ( x `H.outer` y
    , \d -> ( d `H.app` y
            , H.tr d `H.app` x)
    )
{-# INLINE outer #-}

-- | Note: if possible, use the potentially much more performant
-- 'zipWithVector''.
zipWithVector
    :: (KnownNat n, Reifies s W)
    => (BVar s H.ℝ -> BVar s H.ℝ -> BVar s H.ℝ)
    -> BVar s (H.R n)
    -> BVar s (H.R n)
    -> BVar s (H.R n)
zipWithVector f x y = isoVar (H.vecR . SVG.convert . runABP)
                             (ABP . SVG.convert . H.rVec)
                    $ B.liftA2 @(ABP (SV.Vector _)) f (iv x) (iv y)
  where
    iv = isoVar (ABP . SVG.convert . H.rVec) (H.vecR . SVG.convert . runABP)
{-# INLINE zipWithVector #-}

zipWithVector'
    :: ( Num (vec n)
       , Backprop (vec n)
       , Storable field
       , Backprop field
       , H.Sized field (vec n) HU.Vector
       , Reifies s W
       )
    => (forall s'. Reifies s' W => BVar s' field -> BVar s' field -> BVar s' field)
    -> BVar s (vec n)
    -> BVar s (vec n)
    -> BVar s (vec n)
zipWithVector' f = liftOp2 . op2 $ \(VG.convert.H.extract->x) (VG.convert.H.extract->y) ->
    let (z, dx, dy) = V.unzip3 $ V.zipWith (\x' -> retup . backprop2 f x') x y
    in  ( fromJust (H.create (VG.convert z))
        , \d -> ( d * fromJust (H.create (VG.convert dx))
                , d * fromJust (H.create (VG.convert dy))
                )
        )
  where
    retup (x, (y, z)) = (x, y, z)
{-# INLINE zipWithVector' #-}

-- | A version of 'zipWithVector'' that is potentially less performant but
-- is based on 'H.zipWithVector' from 'H.Domain'.
dzipWithVector
    :: ( KnownNat n
       , H.Domain field vec mat
       , Num (vec n)
       , Backprop (vec n)
       , Backprop field
       , Reifies s W
       )
    => (forall s'. Reifies s' W => BVar s' field -> BVar s' field -> BVar s' field)
    -> BVar s (vec n)
    -> BVar s (vec n)
    -> BVar s (vec n)
dzipWithVector f = liftOp2 . op2 $ \x y ->
    ( H.zipWithVector (evalBP2 f) x y
    , \d -> let dx = H.zipWithVector (\x' -> fst . gradBP2 f x') x y
                dy = H.zipWithVector (\x' -> snd . gradBP2 f x') x y
            in  (d * dx, d * dy)
    )
{-# INLINE dzipWithVector #-}

det :: ( KnownNat n
       , Num (mat n n)
       , Backprop (mat n n)
       , H.Domain field vec mat
       , H.Sized field (mat n n) d
       , HU.Transposable (mat n n) (mat n n)
       , Reifies s W
       )
    => BVar s (mat n n)
    -> BVar s field
det = liftOp1 . op1 $ \x ->
    let xDet = H.det x
        xInv = H.inv x
    in  ( xDet, \d -> H.konst (d * xDet) * H.tr xInv )
{-# INLINE det #-}

-- | The inverse and the natural log of the determinant together.  If you
-- know you don't need the inverse, it is best to use 'lndet'.
invlndet
    :: forall n mat field vec d s.
       ( KnownNat n
       , Num (mat n n)
       , H.Domain field vec mat
       , H.Sized field (mat n n) d
       , HU.Transposable (mat n n) (mat n n)
       , Backprop field
       , Backprop (mat n n)
       , Reifies s W
       )
    => BVar s (mat n n)
    -> (BVar s (mat n n), (BVar s field, BVar s field))
invlndet v = (t ^^. _1, (t ^^. _2, t ^^. _3))
  where
    o :: Op '[mat n n] (mat n n, field, field)
    o = op1 $ \x ->
      let (i,(ldet, s)) = H.invlndet x
          iTr           = H.tr i
      in  ( (i, ldet, s)
          , \(dI, dLDet, _) ->
                let gradI    = - iTr `H.mul` dI `H.mul` iTr
                    gradLDet = H.konst dLDet * H.tr i
                in  gradI + gradLDet
          )
    {-# INLINE o #-}
    t = liftOp1 o v
    {-# NOINLINE t #-}
{-# INLINE invlndet #-}

-- | The natural log of the determinant.
lndet
    :: forall n mat field vec d s.
       ( KnownNat n
       , Num (mat n n)
       , Backprop (mat n n)
       , H.Domain field vec mat
       , H.Sized field (mat n n) d
       , HU.Transposable (mat n n) (mat n n)
       , Reifies s W
       )
    => BVar s (mat n n)
    -> BVar s field
lndet = liftOp1 . op1 $ \x ->
          let (i,(ldet,_)) = H.invlndet x
          in  (ldet, (* H.tr i) . H.konst)
{-# INLINE lndet #-}

inv :: ( KnownNat n
       , Num (mat n n)
       , Backprop (mat n n)
       , H.Domain field vec mat
       , HU.Transposable (mat n n) (mat n n)
       , Reifies s W
       )
    => BVar s (mat n n)
    -> BVar s (mat n n)
inv = liftOp1 . op1 $ \x ->
    let xInv   = H.inv x
        xInvTr = H.tr xInv
    in  ( xInv, \d -> - xInvTr `H.mul` d `H.mul` xInvTr )
{-# INLINE inv #-}

toRows
    :: forall m n s. (KnownNat m, KnownNat n, Reifies s W)
    => BVar s (H.L m n)
    -> SV.Vector m (BVar s (H.R n))
toRows = runABP . sequenceVar
       . isoVar (coerce (H.lRows @m)) (coerce H.rowsL)
{-# INLINE toRows #-}

toColumns
    :: forall m n s. (KnownNat m, KnownNat n, Reifies s W)
    => BVar s (H.L m n)
    -> SV.Vector n (BVar s (H.R m))
toColumns = runABP . sequenceVar
          . isoVar (coerce (H.lCols @_ @n)) (coerce H.colsL)
{-# INLINE toColumns #-}

fromRows
    :: forall m n s. (KnownNat m, KnownNat n, Reifies s W)
    => SV.Vector m (BVar s (H.R n))
    -> BVar s (H.L m n)
fromRows = isoVar (coerce H.rowsL) (coerce (H.lRows @m))
         . collectVar . ABP
{-# INLINE fromRows #-}

fromColumns
    :: forall m n s. (KnownNat n, KnownNat m, Reifies s W)
    => SV.Vector n (BVar s (H.R m))
    -> BVar s (H.L m n)
fromColumns = isoVar (coerce H.colsL) (coerce (H.lCols @_ @n))
            . collectVar . ABP
{-# INLINE fromColumns #-}

konst
    :: forall t s d q.
     ( H.Sized t s d
     , HU.Container d t
     , Backprop t
     , Reifies q W
     )
    => BVar q t
    -> BVar q s
konst = liftOp1 . op1 $ \x ->
    ( H.konst x
    , HU.sumElements . H.extract
    )
{-# INLINE konst #-}

sumElements
    :: forall t s d q.
     ( H.Sized t s d
     , HU.Container d t
     , Backprop s
     , Reifies q W
     )
    => BVar q s
    -> BVar q t
sumElements = liftOp1 . op1 $ \x ->
    ( HU.sumElements . H.extract $ x
    , H.konst
    )
{-# INLINE sumElements #-}

-- | If there are extra items in the total derivative, they are dropped.
-- If there are missing items, they are treated as zero.
extractV
    :: forall t s q.
       ( H.Sized t s HU.Vector
       , HU.Konst t Int HU.Vector
       , HU.Container HU.Vector t
       , Backprop s
       , Reifies q W
       )
    => BVar q s
    -> BVar q (HU.Vector t)
extractV = liftOp1 . op1 $ \x ->
    let n = H.size x
    in  ( H.extract x
        , \d -> let m  = HU.size d
                    m' = case compare n m of
                            LT -> HU.subVector 0 n d
                            EQ -> d
                            GT -> HU.vjoin [d, HU.konst 0 (n - m)]
                in  fromJust . H.create $ m'
        )
{-# INLINE extractV #-}

-- | If there are extra items in the total derivative, they are dropped.
-- If there are missing items, they are treated as zero.
extractM
    :: forall t s q.
       ( H.Sized t s HU.Matrix
       , Backprop s
       , HU.Konst t (Int, Int) HU.Matrix
       , HU.Container HU.Matrix t
       , Num (HU.Matrix t)
       , Reifies q W
       )
    => BVar q s
    -> BVar q (HU.Matrix t)
extractM = liftOp1 . op1 $ \x ->
    let (xI,xJ) = H.size x
    in  ( H.extract x
        , \d -> let (dI,dJ) = HU.size d
                    m' = case (compare xI dI, compare xJ dJ) of
                           (LT, LT) -> d HU.?? (HU.Take xI, HU.Take xJ)
                           (LT, EQ) -> d HU.?? (HU.Take xI, HU.All)
                           (LT, GT) -> d HU.?? (HU.Take xI, HU.All)
                                HU.||| HU.konst 0 (xI, xJ - dJ)
                           (EQ, LT) -> d HU.?? (HU.All    , HU.Take xJ)
                           (EQ, EQ) -> d
                           (EQ, GT) -> d HU.?? (HU.All, HU.All)
                                HU.||| HU.konst 0 (xI, xJ - dJ)
                           (GT, LT) -> d HU.?? (HU.All, HU.Take xJ)
                                HU.=== HU.konst 0 (xI - dI, xJ)
                           (GT, EQ) -> d HU.?? (HU.All, HU.All)
                                HU.=== HU.konst 0 (xI - dI, xJ)
                           (GT, GT) -> HU.fromBlocks
                              [[d,0                            ]
                              ,[0,HU.konst 0 (xI - dI, xJ - dJ)]
                              ]
                in  fromJust . H.create $ m'
        )
{-# INLINE extractM #-}

create
    :: (H.Sized t s d, Backprop s, Num (d t), Backprop (d t), Reifies q W)
    => BVar q (d t)
    -> Maybe (BVar q s)
create = sequenceVar . isoVar H.create (maybe 0 H.extract)
{-# INLINE create #-}

takeDiag
    :: ( KnownNat n
       , H.Diag (mat n n) (vec n)
       , H.Domain field vec mat
       , Num field
       , Backprop (mat n n)
       , Reifies s W
       )
    => BVar s (mat n n)
    -> BVar s (vec n)
takeDiag = liftOp1 . op1 $ \x ->
    ( H.takeDiag x
    , H.diagR 0
    )
{-# INLINE takeDiag #-}

-- |
-- \[
-- \frac{1}{2} (M + M^T)
-- \]
sym :: (KnownNat n, Reifies s W)
    => BVar s (H.Sq n)
    -> BVar s (H.Sym n)
sym = liftOp1 . op1 $ \x ->
    ( H.sym x
    , H.unSym . H.sym . H.unSym
    )
{-# INLINE sym #-}

-- |
-- \[
-- M^T M
-- \]
mTm :: (KnownNat m, KnownNat n, Reifies s W)
    => BVar s (H.L m n)
    -> BVar s (H.Sym n)
mTm = liftOp1 . op1 $ \x ->
    ( H.mTm x
    , \d -> 2 * (x H.<> H.unSym d)
    )
{-# INLINE mTm #-}

-- | Warning: the gradient is going necessarily symmetric, and so is /not/
-- meant to be used directly.  Rather, it is meant to be used in the middle
-- (or at the end) of a longer computation.
unSym
    :: (KnownNat n, Reifies s W)
    => BVar s (H.Sym n)
    -> BVar s (H.Sq n)
unSym = isoVar H.unSym unsafeCoerce
{-# INLINE unSym #-}

-- | Unicode synonym for '<.>>'
(<·>)
    :: (KnownNat n, Reifies s W)
    => BVar s (H.R n)
    -> BVar s (H.R n)
    -> BVar s H.ℝ
(<·>) = dot
infixr 8 <·>
{-# INLINE (<·>) #-}

afSV :: Backprop a => BE.AddFunc (SV.Vector n a)
afSV = BE.AF (SV.zipWith add)
{-# INLINE afSV #-}

-- support for complex types
toComplex :: (KnownNat n, Reifies s W)
          => BVar s (H.R n) -> BVar s (H.R n) -> BVar s (H.C n)
toComplex = isoVar2 (curry H.toComplex) H.fromComplex
{-# INLINE toComplex #-}

fromComplex :: (KnownNat n, Reifies s W)
            => BVar s (H.C n) -> (BVar s (H.R n), BVar s (H.R n))
fromComplex c = let ri = isoVar H.fromComplex H.toComplex c in (ri ^^. _1, ri ^^. _2)
{-# INLINE fromComplex #-}

complex :: (KnownNat n, Reifies s W)
       => BVar s (H.R n) -> BVar s (H.C n)
complex = isoVar H.complex H.real
{-# INLINE complex #-}

real :: (KnownNat n, Reifies s W)
     => BVar s (H.C n) -> BVar s (H.R n)
real = isoVar H.real (\r -> H.toComplex (r, H.konst 0))
{-# INLINE real #-}

imag :: (KnownNat n, Reifies s W)
     => BVar s (H.C n) -> BVar s (H.R n)
imag = isoVar H.imag (\r -> H.toComplex (H.konst 0, r))
{-# INLINE imag #-}

sqMagnitude :: (KnownNat n, Reifies s W)
            => BVar s (H.C n) -> BVar s (H.R n)
sqMagnitude c = let (r,i) = fromComplex c in r**2 + i**2
{-# INLINE sqMagnitude #-}

magnitude :: (KnownNat n, Reifies s W)
            => BVar s (H.C n) -> BVar s (H.R n)
magnitude c = sqrt $ sqMagnitude c
{-# INLINE magnitude #-}

max_stars_repo_name: mstksg/hmatrix-backprop

max_issues_repo_name: mstksg/hmatrix-backprop

max_forks_repo_name: mstksg/hmatrix-backprop


text:
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE TemplateHaskell      #-}
{-# LANGUAGE TypeSynonymInstances #-}

module PLA
  ( Weight(..)
  , Point'(..)
  , viewX
  , trainPLA
  , trainPLAWithInitial
  , isMisclassified
  ) where

import           Control.Lens
import           Control.Lens.Combinators
import           Control.Lens.Operators
import           Control.Monad.Reader
import           Data.List
import           Data.Traversable           (for)
import           Numeric.LinearAlgebra.Data
import           System.Random

--------------------------------
--  Solutions to homework 1, problems 7-10 of
--  the 'Learning from Data' course
-- http://work.caltech.edu/homeworks.html
-------------------------------------
-------------------------------
-- Types and Typeclass Instances
--------------------------------
type Weight = (R, R, R)

type Point = (R, R)

-- | To be used at a later stage - currently just a test
data Point' =
  Point'
    { _x :: R
    , _y :: R
    }

makeLenses ''Point'

instance Num Weight where
  (+) (x0, x1, x2) (w0, w1, w2) = (w0 + x0, w1 + x1, w2 + x2)
  (-) (x0, x1, x2) (w0, w1, w2) = (x0 - w0, x1 - w1, x2 - w2)
  (*) (x0, x1, x2) (w0, w1, w2) =
    (w2 * x1 - w1 * x2, w2 * x0 - w0 * x2, w1 * x0 - w0 * x1)
  negate (w0, w1, w2) = (negate w0, negate w1, negate w2)
  abs (w0, w1, w2) = (abs w0, abs w1, abs w2)
  fromInteger x = (fromInteger x, 0.0, 0.0)
  signum (w0, w1, w2) = signum (w0, w1, w2)

--------------------------------------------
-- Creating the line for the target function
--------------------------------------------
m :: Point -> Point -> R
m (x1, x2) (z1, z2) = (x2 - z2) / (x1 - z1)

b :: Point -> Point -> R
b (x1, x2) (z1, z2) = x2 - x1 * (x2 - z2) / (x1 - z1)

-- | Takes 2 Points and returns the function of a line connecting them
line :: Point -> Point -> (R -> R)
line point1 point2 = \x -> m point1 point2 * x + b point1 point2

-----------------------------------
-- functions for the Weight type
-----------------------------------
-- | The dot product
dotProd :: Weight -> Weight -> R
dotProd (x0, x1, x2) (w0, w1, w2) = w0 * x0 + w1 * x1 + w2 * x2

-- | Scalar multiplication of the weight vecotr
prod :: R -> Weight -> Weight
prod lambda (x0, x1, x2) = (x0 * lambda, x1 * lambda, x2 * lambda)

-- | Converting a point to a weight by adding the first coordinate (1)
point2weight :: Point -> Weight
point2weight (x1, x2) = (1, x1, x2)

-- | looks up the label of the point by applying the target function
yn :: (R -> R) -> Point -> R
yn yLine (x1, x2) =
  if yLine x1 < x2
    then (-1)
    else 1

-- | Checking if a point is misclassified
isMisclassified yLine weight point = classified /= label
  where
    classified = signum $ dotProd (point2weight point) weight
    label = yn yLine point

updateWeight' yLine w xn = w + prod (yn yLine xn) (point2weight xn)

-- | Performs 1 training epoch (one run through all training points)
-- | performs the perceptron learning algorithm
pla yLine w niter gen xs
  | misclassifiedXs == [] = (w, niter)
  | otherwise = pla yLine wUpdated (niter + 1) nextRandom xs
  where
    wUpdated = updateWeight' yLine w $ misclassifiedXs !! number
    misclassifiedXs = filter (isMisclassified yLine w) xs
    randomRange = (0, length misclassifiedXs - 1)
    (number, nextRandom) = randomR randomRange gen

-- | stores the target function and performs training, returning the weight and the number of epochs performed before convergence
weightReader' :: [Point] -> StdGen -> Reader (R -> R) (Weight, Int)
weightReader' points gen = do
  targetFunc <- ask
  let weightsEpochs = pla targetFunc (0, 0, 0) 0 gen points
  return weightsEpochs

-- | stores the target function and performs training, returning the weight and the number of epochs performed before convergence
weightReader :: [Point] -> StdGen -> Weight -> Reader (R -> R) (Weight, Int)
weightReader points gen initialWeights = do
  targetFunc <- ask
  let weightsEpochs = pla targetFunc initialWeights 0 gen points
  return weightsEpochs

-- |Creates a random Point between (-1,-1) and (1,1)
createRandomPoint = do
  rand1 :: R <- randomRIO (-1, 1)
  rand2 :: R <- randomRIO (-1, 1)
  return (rand1, rand2)

-- |Creates a list of n random Points between (-1,-1) and (1,1)
createRandomPoints n = do
  randList1 :: [R] <- for [1 .. n] $ \_ -> randomRIO (-1, 1)
  randList2 :: [R] <- for [1 .. n] $ \_ -> randomRIO (-1, 1)
  return $ zip randList1 randList2-- | creates an infinite list of permutations of the list provided as an argument

pointPermutations = cycle . permutations

-- | Just a test for the point lens
viewX point = view x point

---------------------------------------------------------------
-- The final loop that creates the target function and training and testing data,
-- performs training and testing
-------------------------------------------------------------------
-- | performs the PLA
trainPLA :: Int -> IO ()
trainPLA nPoints = do
  point1 <- createRandomPoint
  point2 <- createRandomPoint
  stdgen <- getStdGen
  let target = line point1 point2
  trainpoints <- createRandomPoints nPoints
  let (weights, epochs) = runReader (weightReader' trainpoints stdgen) target
  putStr "Weights: "
  print weights
  putStr "Epochs: "
  print epochs
  testpoints <- createRandomPoints 1000
  let misclassifiedTestPoints =
        length $ filter (isMisclassified target weights) testpoints
  let prob = fromIntegral misclassifiedTestPoints / 1000
  putStr "Probability f /= g: "
  print prob

-- | performs the PLA with a set of initial weights, training points and a target function
trainPLAWithInitial :: (R, R, R) -> [Point] -> (R -> R) -> IO (Weight, Int)
trainPLAWithInitial initialWeights trainPoints target = do
  stdgen <- getStdGen
  let (weights, epochs) =
        runReader (weightReader trainPoints stdgen initialWeights) target
  return (weights, epochs)


text:
module Main
  ( main
  ) where

import Advent (getInputLines)
import Data.Complex (Complex(..))
import Data.List (foldl')

main :: IO ()
main =
  do i <- getInputLines parse 12
     print (part1 i)
     print (part2 i)

type C = Complex Double

parse :: String -> (Char,C)
parse (x:(read -> n)) = (x,n:+0)
parse xs = error (show xs)

dir, rot :: Char -> Maybe C

dir 'N' = Just (   0  :+   1  )
dir 'E' = Just (   1  :+   0  )
dir 'S' = Just (   0  :+ (-1) )
dir 'W' = Just ( (-1) :+   0  )
dir  _  = Nothing

rot 'L' = Just (   0  :+   1  )
rot 'R' = Just (   0  :+ (-1) )
rot  _  = Nothing

part1, part2 :: [(Char,C)] -> Int
step1, step2 :: (C,C) -> (Char,C) -> (C,C)

part1 ins = abs (round x) + abs (round y)
  where
    x :+ y = fst $ foldl' step1 ( 0 :+ 0 , 1 :+ 0 ) ins

step1 (p,h) (a,n)
  | Just v <- dir a = ( p + v*n , h             )
  | Just v <- rot a = ( p       , h * v**(n/90) )
  | otherwise       = ( p + h*n , h             )

part2 ins = abs (round x) + abs (round y)
  where
    x :+ y = fst $ foldl' step2 ( 0 :+ 0 , 10 :+ 1 ) ins

step2 (p,w) (a,n)
  | Just v <- dir a = ( p       , w + v*n       )
  | otherwise       = step1 (p,w) (a,n)


text:
{-# LANGUAGE CApiFFI #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE UndecidableInstances #-}

module LatticeSymmetries.CSR where

-- import Data.Binary (Binary (..))

-- import Data.Vector.Binary

-- import Control.Exception.Safe (bracket, impureThrow, throwIO)
import Control.Monad.Primitive (PrimMonad, PrimState)
import Control.Monad.ST
-- import qualified Control.Monad.ST.Unsafe (unsafeIOToST)
-- import Data.Aeson
-- import Data.Aeson.Types (typeMismatch)

-- import Data.Scientific (toRealFloat)

import Data.Bits (Bits, toIntegralSized)
import Data.Complex
import qualified Data.List
import Data.Type.Equality
-- import qualified Data.List.NonEmpty as NonEmpty

-- import qualified Data.Vector.Fusion.Stream.Monadic as Stream

-- import qualified Data.Vector.Storable as S
import qualified Data.Vector
import qualified Data.Vector as B
import qualified Data.Vector.Algorithms.Intro as Intro
import Data.Vector.Fusion.Bundle (Bundle)
import qualified Data.Vector.Fusion.Bundle as Bundle (inplace)
import qualified Data.Vector.Fusion.Bundle.Monadic as Bundle
import Data.Vector.Fusion.Bundle.Size (Size (..), toMax)
import Data.Vector.Fusion.Stream.Monadic (Step (..), Stream (..))
import qualified Data.Vector.Fusion.Util (unId)
import Data.Vector.Generic ((!))
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Generic.Mutable as GM
import qualified Data.Vector.Storable as S
import qualified Data.Vector.Unboxed
-- import qualified Data.Vector.Storable.Mutable as SM
-- import Data.Vector.Unboxed (Unbox)
-- import qualified Data.Vector.Unboxed as U
-- import Data.Yaml (decodeFileWithWarnings)
-- import Foreign.C.String (CString, peekCString)
import Foreign.C.Types (CInt (..), CUInt (..))
-- import Foreign.ForeignPtr
-- import Foreign.ForeignPtr.Unsafe (unsafeForeignPtrToPtr)
-- import Foreign.Marshal.Alloc (alloca, free, malloc, mallocBytes)
-- import Foreign.Marshal.Array (newArray, withArrayLen)
import Foreign.Marshal.Utils (with)
import Foreign.Ptr (FunPtr, Ptr, castPtr)
-- import Foreign.StablePtr
import Foreign.Storable (Storable (..))
import qualified GHC.Exts as GHC (IsList (..))
import GHC.TypeLits
import qualified GHC.TypeLits as GHC
-- import qualified GHC.ForeignPtr as GHC (Finalizers (..), ForeignPtr (..), ForeignPtrContents (..))
-- import GHC.Generics
-- import qualified GHC.IORef as GHC (atomicSwapIORef)
-- import GHC.Prim
-- import qualified GHC.Ptr as GHC (Ptr (..))
-- import qualified Language.C.Inline as C
-- import qualified Language.C.Inline.Unsafe as CU
import LatticeSymmetries.Dense
-- import LatticeSymmetries.IO
-- import LatticeSymmetries.Types
import qualified System.IO.Unsafe
import qualified Unsafe.Coerce
-- import qualified System.Mem.Weak
import Prelude hiding (group, product, sort)

-- typedef struct ls_csr_matrix {
--   unsigned *offsets;
--   unsigned *columns;
--   _Complex double *off_diag_elements;
--   _Complex double *diag_elements;
--   unsigned dimension;
--   unsigned number_nonzero;
-- } ls_csr_matrix;
data Ccsr_matrix = Ccsr_matrix
  { c_csr_matrix_offsets :: !(Ptr CUInt),
    c_csr_matrix_columns :: !(Ptr CUInt),
    c_csr_matrix_off_diag_elements :: !(Ptr (Complex Double)),
    c_csr_matrix_diag_elements :: !(Ptr (Complex Double)),
    c_csr_matrix_dimension :: !CUInt,
    c_csr_matrix_number_nonzero :: !CUInt
  }

instance Storable Ccsr_matrix where
  alignment _ = 8
  sizeOf _ = 40
  peek p =
    Ccsr_matrix
      <$> peekByteOff p 0
      <*> peekByteOff p 8
      <*> peekByteOff p 16
      <*> peekByteOff p 24
      <*> peekByteOff p 32
      <*> peekByteOff p 36
  poke p x = do
    pokeByteOff p 0 (c_csr_matrix_offsets x)
    pokeByteOff p 8 (c_csr_matrix_columns x)
    pokeByteOff p 16 (c_csr_matrix_off_diag_elements x)
    pokeByteOff p 24 (c_csr_matrix_diag_elements x)
    pokeByteOff p 32 (c_csr_matrix_dimension x)
    pokeByteOff p 36 (c_csr_matrix_number_nonzero x)

data CsrMatrix
  = CsrMatrix
      !(S.Vector CUInt)
      !(S.Vector CUInt)
      !(S.Vector (Complex Double))
      !(S.Vector (Complex Double))
  deriving stock (Show, Eq)

withCsrMatrix :: CsrMatrix -> (Ptr Ccsr_matrix -> IO a) -> IO a
withCsrMatrix m@(CsrMatrix offsets columns offDiagElems diagElems) action =
  S.unsafeWith offsets $ \offsetsPtr ->
    S.unsafeWith columns $ \columnsPtr ->
      S.unsafeWith offDiagElems $ \offDiagElemsPtr ->
        S.unsafeWith diagElems $ \diagElemsPtr ->
          with (Ccsr_matrix offsetsPtr columnsPtr offDiagElemsPtr diagElemsPtr c_dimension c_number_non_zero) action
  where
    c_dimension = fromIntegral (csrDim m)
    c_number_non_zero = fromIntegral (csrNumberNonZero m)

unsafeNewCsrMatrix :: Int -> Int -> IO CsrMatrix
unsafeNewCsrMatrix dimension nnz =
  pure $
    CsrMatrix
      (G.replicate (dimension + 1) 0)
      (G.replicate nnz 0)
      (G.replicate nnz 0)
      (G.replicate dimension 0)

csrMatrixShrink :: HasCallStack => Int -> CsrMatrix -> CsrMatrix
csrMatrixShrink nnz m@(CsrMatrix offsets columns offDiagElems diagElems)
  | nnz == csrNumberNonZero m = m
  | nnz < csrNumberNonZero m = CsrMatrix offsets (G.take nnz columns) (G.take nnz offDiagElems) diagElems
  | otherwise = error "general resizing not implemented"

csrDim :: CsrMatrix -> Int
csrDim (CsrMatrix _ _ _ diagElems) = G.length diagElems

csrNumberNonZero :: CsrMatrix -> Int
csrNumberNonZero (CsrMatrix _ _ offDiagElems _) = G.length offDiagElems

csrMatrixFromDense ::
  (HasCallStack, G.Vector v (Complex Double), G.Vector v Int) =>
  DenseMatrix v (Complex Double) ->
  CsrMatrix
csrMatrixFromDense m@(DenseMatrix nRows nCols v)
  | nRows == nCols = runST $ do
    let nnz = countOffDiagNonZero m
        diagElems = G.convert $ extractDiagonal m
    offsetsBuffer <- GM.new (nRows + 1)
    columnsBuffer <- GM.new nnz
    offDiagElemsBuffer <- GM.new nnz
    let go2 !i !j !offset
          | i == j = go2 i (j + 1) offset
          | j < nCols =
            let x = indexDenseMatrix m (i, j)
             in if x /= 0
                  then do
                    GM.write columnsBuffer offset (fromIntegral j)
                    GM.write offDiagElemsBuffer offset x
                    go2 i (j + 1) (offset + 1)
                  else go2 i (j + 1) offset
          | otherwise = pure offset
        go1 !i !offset
          | i < nRows = do
            offset' <- go2 i 0 offset
            GM.write offsetsBuffer (i + 1) (fromIntegral offset')
            go1 (i + 1) offset'
          | otherwise = pure offset
    GM.write offsetsBuffer 0 0
    offset <- go1 0 0
    unless (nnz == offset) $ error "this is probably a bug"
    CsrMatrix
      <$> G.unsafeFreeze offsetsBuffer
      <*> G.unsafeFreeze columnsBuffer
      <*> G.unsafeFreeze offDiagElemsBuffer
      <*> pure diagElems
  | otherwise = error "expected a square matrix"

foreign import capi unsafe "csr.h ls_csr_plus"
  ls_csr_plus :: Ptr Ccsr_matrix -> Ptr Ccsr_matrix -> Ptr Ccsr_matrix -> IO ()

foreign import capi unsafe "csr.h ls_csr_minus"
  ls_csr_minus :: Ptr Ccsr_matrix -> Ptr Ccsr_matrix -> Ptr Ccsr_matrix -> IO ()

foreign import capi unsafe "csr.h ls_csr_times"
  ls_csr_times :: Ptr Ccsr_matrix -> Ptr Ccsr_matrix -> Ptr Ccsr_matrix -> IO ()

foreign import capi unsafe "csr.h ls_csr_kron"
  ls_csr_kron :: Ptr Ccsr_matrix -> Ptr Ccsr_matrix -> Ptr Ccsr_matrix -> IO ()

binaryOp :: (Ptr Ccsr_matrix -> Ptr Ccsr_matrix -> Ptr Ccsr_matrix -> IO ()) -> CsrMatrix -> CsrMatrix -> CsrMatrix
binaryOp opKernel a b = System.IO.Unsafe.unsafePerformIO $ do
  let pessimisticNNZ = csrNumberNonZero a + csrNumberNonZero b
  c <- unsafeNewCsrMatrix (csrDim a) pessimisticNNZ
  nnz <-
    withCsrMatrix a $ \aPtr ->
      withCsrMatrix b $ \bPtr ->
        withCsrMatrix c $ \cPtr -> do
          opKernel aPtr bPtr cPtr
          fromIntegral . c_csr_matrix_number_nonzero <$> peek cPtr
  pure $ csrMatrixShrink nnz c

instance Num CsrMatrix where
  (+) = binaryOp ls_csr_plus
  (-) = binaryOp ls_csr_minus
  (*) = binaryOp ls_csr_times
  abs (CsrMatrix offsets columns offDiagElems diagElems) = CsrMatrix offsets columns (G.map abs offDiagElems) (G.map abs diagElems)
  signum (CsrMatrix offsets columns offDiagElems diagElems) = CsrMatrix offsets columns (G.map signum offDiagElems) (G.map signum diagElems)
  fromInteger _ = error "Num instance of CsrMatrix does not implement fromInteger"

csrScale :: Complex Double -> CsrMatrix -> CsrMatrix
csrScale c (CsrMatrix offsets columns offDiagElems diagElems) =
  CsrMatrix offsets columns (G.map (c *) offDiagElems) (G.map (c *) diagElems)

csrKron :: CsrMatrix -> CsrMatrix -> CsrMatrix
csrKron a b = System.IO.Unsafe.unsafePerformIO $ do
  let dimension = csrDim a * csrDim b
      nnz = csrNumberNonZero a * csrNumberNonZero b
  c <- unsafeNewCsrMatrix dimension nnz
  nnz' <-
    withCsrMatrix a $ \aPtr ->
      withCsrMatrix b $ \bPtr ->
        withCsrMatrix c $ \cPtr -> do
          ls_csr_kron aPtr bPtr cPtr
          fromIntegral . c_csr_matrix_number_nonzero <$> peek cPtr
  unless (nnz == nnz') $
    error $ "this is probably a bug: " <> show nnz <> " /= " <> show nnz'
  pure c

csrKronMany :: HasCallStack => [CsrMatrix] -> CsrMatrix
csrKronMany [] = error "expected a non-empty list of matrices"
csrKronMany xs = Data.List.foldl1' csrKron xs

max_stars_repo_name: twesterhout/lattice-symmetries-haskell

max_issues_repo_name: twesterhout/lattice-symmetries-haskell

max_forks_repo_name: twesterhout/lattice-symmetries-haskell


text:
{-# LANGUAGE ScopedTypeVariables, TypeFamilies #-}

module Fuml.Base.KNN where

import qualified Data.Vector.Storable as VS
import Numeric.LinearAlgebra
import Data.List (nub, sortBy)
import Data.Ord (comparing)

euclideanDistance :: Vector Double -> Vector Double -> Double
euclideanDistance v1 v2 = sqrt $ VS.sum $ VS.map (^2) $ VS.zipWith (-) v1 v2

weightedBoolVote :: [(Double, Bool)] -> Double
weightedBoolVote distBools =
  let wtrue = sum $ map (recip . fst) $ filter snd distBools
      wfalse = sum $ map (recip . fst) $ filter (not . snd) distBools
  in exp wtrue / (exp wtrue + exp wfalse)

majorityVote :: Eq a => [(Double, a)] -> a
majorityVote distXs = let classes = nub $ map snd distXs
                          occurences c = (c,negate $ length $ filter ((==c) . snd) distXs)
                      in fst $ head $ sortBy (comparing snd) $ map occurences classes

weightedMajorityVote :: Eq a => [(Double, a)] -> a
weightedMajorityVote distXs
  = let classes = nub $ map snd distXs
        weight c = (c,negate $ sum $ map (recip . fst) $ filter ((==c) . snd) distXs)
        in fst $ head $ sortBy (comparing snd) $ map weight classes

avgVote :: [(Double, Double)] -> Double
avgVote distXs = let n = realToFrac $ length distXs
                 in (sum $ map snd distXs) / n

weightedAvgVote :: [(Double, Double)] -> Double
weightedAvgVote distXs
  = let wsum = sum $ map (recip . fst) distXs
    in  (sum $ map (uncurry (*)) distXs) / wsum

max_stars_repo_name: ekalosak/open

max_issues_repo_name: ekalosak/open

max_forks_repo_name: ekalosak/open


text:
{- -}

module Neural.Digits where

import Control.DeepSeq
import Control.Monad
import Numeric.LinearAlgebra
import System.Random

import Neural.FF
import Neural.LayerSpawners
import Neural.MNIST

-- 1st layer fn is a dummy. I might replace with a dedicated dummy layer fn
spawnDigitClassifierNet :: IO (FFNetwork)
spawnDigitClassifierNet = spawnFFNetwork [(fullyConnectedLayer, (0,0), (28, 28), layerOpts),
                                          (convLayer, (28, 28), (5,5), layerOpts {
                                              layerDepth = 5,
                                              pooling = maxPooling,
                                              poolingBP = maxPoolingBP,
                                              poolSize = (2,2)
                                              }),
                                          (fullyConnectedLayer, (12, 12), (10, 1), layerOpts)]

trainDigitClassifierNet :: FFNetwork -> [((Matrix Double), Int)] -> [((Matrix Double), Int)] -> Int -> IO (FFNetwork, [(Matrix Double, Int)])
trainDigitClassifierNet net trainingSet testSet nEpochs = go net trainingSet 1
  where go net ts i | i <= nEpochs = do
                        (newNet, newTS) <- getStdRandom (\gen -> sgd net gen ts 9 16 0.1 5.0)
                        let nCorrect = (evaluate newNet testSet 5.0)

                        putStrLn $ "Epoch " ++ (show i) ++ " of " ++ (show nEpochs) ++ ": "
                          ++ (show nCorrect) ++ "/" ++ (show $ length testSet)

                        {-
                        -- print out features
                        mapM_ (\layer ->
                                if (length $ weights layer) > 1
                                then (printMNISTImage $
                                      (\(w_1:w_) ->
                                        foldl (\w w_ -> w ||| (konst 0.0 (5,3) :: Matrix Double) ||| w_) w_1 w_)
                                      $ weights layer) >> putStrLn ""
                                else putStrLn "") $ layers net -}

                        go newNet newTS $ i + 1

                    | otherwise = return (net, trainingSet)

testDigitClassifierNet :: FFNetwork -> [((Matrix Double), Int)] -> IO ()
testDigitClassifierNet net testSet = do
  let tests = take 3 testSet
  putStrLn "Example tests: "
  printMNISTImage (foldl (\mat (input,expected) -> mat ||| input) ((\(x,y) -> x) (head tests)) (tail tests))

  putStrLn $ "Guesses: " ++
    (show $ map (\(input,expected) -> convertOutput $ ffOut net [input]) tests)

  putStrLn $ "Full evaluation: " ++ (show $ evaluate net testSet 1.0) ++ "/" ++ (show $ length testSet)