text:
-- calculates nash equilibriums, does related backend stuff
{-# LANGUAGE OverloadedStrings #-}

module Game
    ( gameSimplePartial, gameSimplePartialOpts, gameSimple
    , calcPureEVs, calcEVSimpleCore, calcVarSimpleCore, calcSDSimpleCore
    , gameComplex, resultComplex
    , solvePure, solvePureBy, solvePureOn, solvePureZeroSum
    , solvePureOptimal, solvePureOptimalZeroSum
    , solveSimpleCoreNaive
    , calcEVSimpleCoreNaive
    , solveSimple, solveComplex
    , resultComplex
    , Game (Game, gameCName, attCName, defCName, gameData, outcomesC)
    , Result (Result, evc, evr, sdc, sdr, weightsAtts, weightsDefs)
    , ResultSimple  (ResultSimple, evSimpleCol, evSimpleRow, sdSimpleCol, sdSimpleRow, weightsColsSimple, weightsRowsSimple)
    ) where

import GameData

import Data.List
import Data.List.Split
import Data.Function
import Data.Maybe
import Data.Text (Text, unpack)
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data
import Data.Map (Map)
import qualified Data.Map as Map
import Safe

type EV = Double
type Var = Double
type SD = Double
type Weights = [Double]
type Coord = (Int,Int)
type RemovedElements = [Int]
type WeightedOption = (Opt, Opt, Double)
type BiWeightedOption = (Opt, Opt, Double, Double)
type PayoffsRaw = [[Double]]
type Payoffs = Matrix Double

gameSimplePartial :: PayoffsRaw -> Payoffs
gameSimplePartial = fromLists

gameSimplePartialOpts :: [WeightedOption] -> PayoffsRaw
gameSimplePartialOpts opts =
    let
        fst3 (x,_,_) = x
        snd3 (_,x,_) = x
        thd3 (_,_,x) = x

        cols = nub . map (fst . fst3) $ opts
        rows = nub . map (fst . snd3) $ opts
        grid = map (map thd3) . transpose . chunksOf (length rows) $ opts
    in
        grid

gameSimple :: Text -> [Text] -> [Text] -> PayoffsRaw -> PayoffsRaw -> GameSimple
gameSimple n c r m1 m2 = GameSimple n c r (fromLists m1) (fromLists m2) $ solveSimple m1 m2

calcPureEVs :: Payoffs -> Weights -> [EV]
calcPureEVs m cols = concat . toLists $ m Numeric.LinearAlgebra.<> (col (fmap (/ sum cols) cols))

calcEVSimpleCore :: Payoffs -> Weights -> Weights -> EV
calcEVSimpleCore m cols rows = (head . head . toLists $ (row rows) Numeric.LinearAlgebra.<> m Numeric.LinearAlgebra.<> (col cols))

calcVarSimpleCore :: Payoffs -> Weights -> Weights -> Var
calcVarSimpleCore mat cols rows = calcEVSimpleCore (cmap (\x -> (x - calcEVSimpleCore mat cols rows)**2) mat) cols rows

calcSDSimpleCore :: Payoffs -> Weights -> Weights -> SD
calcSDSimpleCore a b = sqrt . calcVarSimpleCore a b

gameComplex :: Text -> Text -> Text -> [BiWeightedOption] -> Game
gameComplex t1 t2 t3 m = Game t1 t2 t3 (sort m) . resultComplex m $ solveComplex m

resultComplex :: [(Opt, Opt, Double, Double)] -> ResultSimple -> Result
resultComplex opts (ResultSimple evc evr sdc sdr cw rw) = do
    let colnames = nub . map (fst . fst4) $ opts
    let rownames = nub . map (fst . snd4) $ opts
    Result evc evr sdc sdr (zip colnames cw) (zip rownames rw)

-- from here on out, stuff for solving games

-- the fundamental algorithm, given some ordering method to know which values are better
solvePureBy :: (Eq a) => (a -> a -> Ordering) -> [[a]] -> [[a]] -> [Coord]
solvePureBy f mc mr =
    let
        width = length . head $ mc
        height = length mc

        -- column player: picks the highest columns for each row
        mcIsBestResponse = map (\col -> map (== maximumBy f col) col) mc
        -- row player: picks the highest rows for each coloumn
        mrIsBestResponse = transpose . map (\row -> map (== maximumBy f row) row) . transpose $ mr
        bestResponses = map (uncurry zip) $ zip mcIsBestResponse mrIsBestResponse

        coords = map (\r -> map (\c -> (c,r)) [0..width-1]) [0..height-1]
        zipped = map (uncurry zip) $ zip coords bestResponses
        zippedFiltered = map fst . filter (\((x,y),(c,r)) -> and[c,r]) . concat $ zipped
    in
        zippedFiltered

-- solvePureBy applied to mapped values
solvePureOn :: (Ord c) => (a -> c) -> (b -> c) -> [[a]] -> [[b]] -> [Coord]
solvePureOn f g mc mr = solvePureBy compare (map (map f) mc) (map (map g) mr)

solvePure :: (Ord a) => [[a]] -> [[a]] -> [Coord]
solvePure mc mr = solvePureOn id id mc mr

solvePureZeroSum :: (Ord a, Num a) => [[a]] -> [Coord]
solvePureZeroSum m = solvePureOn id negate m m

evPure :: [[a]] -> Coord -> a
evPure m (c,r) = (m!!c)!!r

solvePureOptimalBy :: (Ord a) => (a -> a -> Ordering) -> ([[a]] -> Coord -> a) -> [[a]] -> [[a]] -> Maybe Coord
solvePureOptimalBy f ev mc mr =
    let
        pures = solvePureBy f mc mr
        -- note that we have already optimised the row choice per column choice, as these are all, by definition, nash equilibria
        -- thus, simply choose the optimal nash equilibrium for column ev
        bestPure = maximumByMay (compare `on` (\p -> ev mc p)) pures
    in
        bestPure

solvePureOptimalOn :: (Ord c) => (a -> c) -> (b -> c) -> ([[c]] -> Coord -> c) -> [[a]] -> [[b]] -> Maybe Coord
solvePureOptimalOn f g ev mc mr = solvePureOptimalBy compare ev (map (map f) mc) (map (map g) mr)

solvePureOptimal :: (Ord a) => [[a]] -> [[a]] -> Maybe Coord
solvePureOptimal mc mr = solvePureOptimalOn id id evPure mc mr

solvePureOptimalZeroSum :: (Ord a, Num a) => [[a]] -> Maybe Coord
solvePureOptimalZeroSum m = solvePureOptimal m . map (map negate) $ m

solveSimpleCoreNaive :: Matrix Double -> Maybe Weights
solveSimpleCoreNaive m =
    do
        case (rows m, cols m) of
            -- not really much choice
            (_,1) -> return [1]
            -- they only have one option, so zip your options to their position, get the best option's position, and pad a [1] with zeroes as appropriate to the position
            (1,c) -> return . (\n -> ((take n (repeat 0)) ++ [1] ++ (take (c - n - 1) (repeat 0)))) . fst . maximumBy (compare `on` snd) . zip [0..] . concat . toLists $ m
            -- use linear algebra to find the weights such that unexploitable play is achieved, ie, the ev of each of the opponent's options is equal
            _ -> let
                        weights = normalise . concat . toLists $ ((pinv m) Numeric.LinearAlgebra.<> (((rows m)><1) (repeat (fromInteger 1))))
                    in
                        if all (>=0) weights
                        then return weights
                        else Nothing
    where
        normalise :: (Fractional a) => [a] -> [a]
        normalise xs = map (/ (sum xs)) xs

calcEVSimpleCoreNaive :: Matrix Double -> Maybe EV
calcEVSimpleCoreNaive m =
    do
        wc <- solveSimpleCoreNaive m
        wr <- solveSimpleCoreNaive . tr' $ m
        return $ calcEVSimpleCore m wc wr

gameSupports :: (Ord a, Num a) => [[a]] -> [[a]] -> [[((RemovedElements, RemovedElements), ([[a]], [[a]]))]]
gameSupports m1 m2 =
    let
        widther = length . head
        heighter = length
        w = widther m1 - 1
        h = heighter m1 - 1

        removeIndexes :: RemovedElements -> [a] -> [a]
        removeIndexes [] xs = xs
        removeIndexes _ [] = error "Can't remove from an empty list!"
        removeIndexes [0] (x:xs) = xs
        removeIndexes [n] (x:xs) = x:(removeIndexes [n-1] xs)
        removeIndexes (n:ns) xs = removeIndexes (fmap (subtract 1) ns) (removeIndexes [n] xs)

        subm :: (RemovedElements,RemovedElements) -> [[a]] -> [[a]]
        subm coords = map (removeIndexes (fst coords)) . removeIndexes (snd coords)

        supportsColsRaw = zip [0..] . init . subsequences $ [0..w] -- [[],[0],[1],...,[n],[0,1],[0,2],...,[0,n],[0,1,2],...]
        supportsRowsRaw = zip [0..] . init . subsequences $ [0..h]
        supportsPairs = map (\(x,c) -> map (\(y,r) -> ((c,r), (subm (c,r) m1, subm (c,r) m2))) supportsRowsRaw) supportsColsRaw
    in
        supportsPairs

gameSupportsZeroSum :: (Ord a, Num a) => [[a]] -> [[((RemovedElements, RemovedElements), ([[a]], [[a]]))]]
gameSupportsZeroSum m = gameSupports m (map (map negate) m)

supportIntoValidatedWeightsEVsSDs :: ((RemovedElements, RemovedElements), ([[Double]], [[Double]])) -> Maybe ((Weights, Weights), ((EV, EV), (SD, SD)))
supportIntoValidatedWeightsEVsSDs ((sc,sr), (mc,mr)) =
    do
        let (g1,g2) = (gameSimplePartial mc, tr' $ gameSimplePartial mr)
        w1 <- solveSimpleCoreNaive g1
        w2 <- solveSimpleCoreNaive g2
        let ws = (addIndexes 0 sc w1, addIndexes 0 sr w2)
        let evs = (calcEVSimpleCore g1 w1 w2, calcEVSimpleCore g2 w2 w1)
        let sds = (calcSDSimpleCore g1 w1 w2, calcSDSimpleCore g2 w2 w1)
        return (ws, (evs, sds))
    where
        addIndexes :: a -> [Int] -> [a] -> [a]
        addIndexes _ [] xs = xs
        addIndexes def [0] xs = def:xs
        addIndexes _ _ [] = error "Add index out of bounds"
        addIndexes def [n] (x:xs) = x:(addIndexes def [n-1] xs)
        addIndexes def (n:ns) xs = addIndexes def ns (addIndexes def [n] xs)

supportsIntoValidatedWeightsEVsSDs :: [[((RemovedElements, RemovedElements), ([[Double]], [[Double]]))]] -> [[Maybe ((Weights, Weights), ((EV, EV), (SD, SD)))]]
supportsIntoValidatedWeightsEVsSDs = map (map supportIntoValidatedWeightsEVsSDs)

solveSimple :: [[Double]] -> [[Double]] -> ResultSimple
solveSimple mc mr =
    let
        sups = gameSupports mc mr
        strats = supportsIntoValidatedWeightsEVsSDs sups

        maximumOn f = maximumBy (compare `on` f)
        colEv = fmap (fst . fst . snd)
        rowEv = fmap (snd . fst . snd)

        stratsBestResponses = map (maximumOn rowEv) strats
        stratsBest = maximumOn colEv stratsBestResponses

        res (Just ((wc,wr),((evc,evr),(sdc,sdr)))) = ResultSimple evc evr sdc sdr wc wr
    in
        res stratsBest

solveComplex :: Ord a => [((a, Maybe Double), (a, Maybe Double), Double, Double)] -> ResultSimple
solveComplex gdata =
    let
        typeCheck ((_,c), (_,r), _, _) = (fromEnum . isNothing $ c) + ((2*) . fromEnum . isNothing $ r)

        pureStrategy :: Int -> Int -> [Double]
        pureStrategy l x = take l ((take (x-1) . repeat $ 0) ++ (1:[0..]))
        --  get all the subgame data and split it according to type: neither = Nothing, col = Nothing, row = Nothing, both = Nothing
        --  evs :: [[((Text, Maybe Double), (Text, Maybe Double), Double)]]
        evs = map (\i -> filter (\outcome -> typeCheck outcome == i) $ gdata) [0..3]
    in
        case evs of
            [both, [], [], []]  -> do -- both are fixed, so we just copy the weights into the ResultSimple
                let rows = nub . map (fst . fst4) $ both
                let cols = nub . map (fst . snd4) $ both
                let grid = transpose . chunksOf (length rows) $ both
                let outsc = map (map thd4) grid
                let outsr = map (map fth4) grid

                let gc = fromLists outsc
                let gr = fromLists outsr
                let cw = map (maybe (error "???") id . snd) . nub . map fst4 $ both
                let rw = map (maybe (error "???") id . snd) . nub . map snd4 $ both
                ResultSimple (calcEVSimpleCore gc cw rw) (calcEVSimpleCore gr cw rw) (calcSDSimpleCore gc cw rw) (calcSDSimpleCore gr cw rw) cw rw

            [[], c, [], []]      -> do -- typeCheck 1 -> columns unfixed, rows fixed, evaluate columns
                let cols = nub . map (fst . fst4) $ c -- each column option, in the form of name::Text
                let rows = map (\(a, Just b) -> (a,b)) . nub . map snd4 $ c -- each row option, in the form of (name, weight)::(Text, Double)
                let grid = transpose . chunksOf (length rows) $ c
                let outsc = map (map thd4) grid
                let outsr = map (map fth4) grid
                let gc = fromLists outsc
                let gr = fromLists outsr

                let colStrats = map (pureStrategy . length $ cols) [0..length cols] -- every pure strategy to try against the fixed row strategy
                let evs = map (\strat -> calcEVSimpleCore gc strat (map snd rows)) colStrats -- every ev

                let (cw, ev) = maximumBy (compare `on` snd) . zip colStrats $ evs -- get the optimal pure strategy for the column player
                let rw = map snd rows
                ResultSimple ev (calcEVSimpleCore gr cw rw) (calcSDSimpleCore gc cw rw) (calcSDSimpleCore gr cw rw) cw rw

            [[], [], r, []]      -> do -- the same as above, but cols fixed rows unfixed
                let cols = map (\(a, Just b) -> (a,b)) . nub . map fst4 $ r
                let rows = nub . map (fst . snd4) $ r
                let grid = transpose . chunksOf (length rows) $ r
                let outsc = map (map thd4) grid
                let outsr = map (map fth4) grid
                let gc = fromLists outsc
                let gr = fromLists outsr

                let rowStrats = map (pureStrategy . length $ rows) [0..length rows]
                let evs = map (\strat -> calcEVSimpleCore gr (map snd cols) strat) rowStrats

                let cw = map snd cols
                let (rw, ev) = minimumBy (compare `on` snd) . zip rowStrats $ evs -- the defender wants to minimise the ev
                ResultSimple (calcEVSimpleCore gc cw rw) ev (calcSDSimpleCore gc cw rw) (calcSDSimpleCore gr cw rw) cw rw

            [[], [], [], neither]     -> do -- both are unfixed, so we just calculate it as a normal GameSimple
                let cols = nub . map (fst . fst4) $ neither
                let rows = nub . map (fst . snd4) $ neither
                let grid = transpose . chunksOf (length rows) $ neither
                let outsc = map (map thd4) grid
                let outsr = map (map fth4) grid

                solveSimple outsc outsr

            _                      -> error "both fixed and unfixed options for a single player" -- no. don't. this makes no sense!!! the only reason i can even think you'd want to try is equivalent to just nesting in another mixup, so just do that.

fst4 (x,_,_,_) = x
snd4 (_,x,_,_) = x
thd4 (_,_,x,_) = x
fth4 (_,_,_,x) = x

max_stars_repo_name: StaccatoSemibreve/Mixup-Analyser

max_issues_repo_name: StaccatoSemibreve/Mixup-Analyser

max_forks_repo_name: StaccatoSemibreve/Mixup-Analyser


text:
{-# LANGUAGE BangPatterns      #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TupleSections     #-}
module Releaser.SettingsRouting
    ( routing
    , bnNbn
    , mapProductType
    , productTypes
    , ConfigRouting (..)
    , ConfigProcTimes (..)
    , queues
    , isRoutedOver'
    , machines
    , allBlocks
    , procTimes
    , procTimesConst
    ) where

import           Control.Arrow                       (second)
import           Data.List                           (nub, sort)
import           Data.Text                           (Text)
import           Prelude
import           SimSim                              hiding (allBlocks, productTypes,
                                                      queues)
import           Statistics.Distribution
import           Statistics.Distribution.Exponential
import           Statistics.Distribution.Uniform
import           System.Random.MWC

import           Releaser.Routing.Ops
import           Releaser.Routing.Type

routing :: ConfigRouting
routing =
  -- routingSingleStageSingleProduct
  -- routingSingleStage
  -- routingDiverging2Stages
  routingDiverging3Stages

procTimes :: ConfigProcTimes
procTimes =
  -- procTimesDiverging2Stages
  -- procTimesDiverging3StagesUnif
  procTimesDiverging3StagesExp

productTypes :: [ProductType]
productTypes = sort $ nub $ map (fst . fst) (configRoutingRoutes routing)

isRoutedOver' :: Block -> ProductType -> Bool
isRoutedOver' bl pt = isRoutedOver (configRoutingRoutes routing) pt bl

bnNbn :: Bool
bnNbn = False -- True

mapProductType :: ProductType -> ProductType
mapProductType x | not bnNbn = x
mapProductType (Product 1) = Product 1
mapProductType (Product 2) = Product 2
mapProductType (Product 3) = Product 1
mapProductType (Product 4) = Product 1
mapProductType (Product 5) = Product 2
mapProductType (Product 6) = Product 1
mapProductType _ = error "function mapProductType not complete in SettingsRouting"

allBlocks :: [Block]
allBlocks = nub $ map (snd . fst) (configRoutingRoutes routing)

queues :: [Block]
queues = filter isQueue allBlocks

machines :: [Block]
machines = filter isMachine allBlocks

--

data ConfigProcTimes = ConfigProcTimes
  { configProcTimesName :: !Text
  , configProcTimes     :: !ProcTimes
  }

test :: IO ()
test = do
  let mkTime = timeFromDouble <$> withSystemRandom (asGenIO $ genContVar (exponential (960 / 80)))
      nr = 10000
  ts <- mapM (const mkTime ) [1..nr]

  print (sum ts / fromIntegral nr)

procTimesDiverging3StagesExp :: ConfigProcTimes
procTimesDiverging3StagesExp = ConfigProcTimes "Proc Times diverging 3 stages: Exp" $
  map (second $ filter ((`elem` productTypes) . fst)) $ filter ((`elem` allBlocks) . fst)
  [ (Machine 1, forAllProducts (fmap timeFromDouble . genContVar (exponential (960  / 80))))
  , (Machine 2, forAllProducts (fmap timeFromDouble . genContVar (exponential (960 / 160))))
  , (Machine 3, forAllProducts (fmap timeFromDouble . genContVar (exponential (960 / 155))))
  , (Machine 4, forAllProducts (fmap timeFromDouble . genContVar (exponential (960 / 210))))
  , (Machine 5, forAllProducts (fmap timeFromDouble . genContVar (exponential (960 / 285))))
  , (Machine 6, forAllProducts (fmap timeFromDouble . genContVar (exponential (960 / 215))))
  ]

  where forAllProducts x = map (,x) productTypes

procTimesDiverging3StagesUnif :: ConfigProcTimes
procTimesDiverging3StagesUnif = ConfigProcTimes "Proc Times diverging 3 stages: Unif" $
  map (second $ filter ((`elem` productTypes) . fst)) $ filter ((`elem` allBlocks) . fst)
  [ (Machine 1, forAllProducts (fmap timeFromDouble . genContVar (uniformDistr (30 / 960) (130 / 960))))
  , (Machine 2, forAllProducts (fmap timeFromDouble . genContVar (uniformDistr (80 / 960) (240 / 960))))
  , (Machine 3, forAllProducts (fmap timeFromDouble . genContVar (uniformDistr (50 / 960) (260 / 960))))
  , (Machine 4, forAllProducts (fmap timeFromDouble . genContVar (uniformDistr (50 / 960) (370 / 960))))
  , (Machine 5, forAllProducts (fmap timeFromDouble . genContVar (uniformDistr (200 / 960) (370 / 960))))
  , (Machine 6, forAllProducts (fmap timeFromDouble . genContVar (uniformDistr (110 / 960) (320 / 960))))
  ]

  where forAllProducts x = map (,x) productTypes

procTimesDiverging2Stages :: ConfigProcTimes
procTimesDiverging2Stages = ConfigProcTimes "Proc Times diverging 2 stages: Unif" $
  map (second $ filter ((`elem` productTypes) . fst)) $ filter ((`elem` allBlocks) . fst)
  [ (Machine 1, [ (Product 1, fmap timeFromDouble . genContVar (uniformDistr (70 / 960) (130 / 960)))
                , (Product 2, fmap timeFromDouble . genContVar (uniformDistr (70 / 960) (130 / 960)))])
  , (Machine 2, [(Product 1, fmap timeFromDouble . genContVar (uniformDistr (130 / 960) (170 / 960)))])
  , (Machine 3, [(Product 2, fmap timeFromDouble . genContVar (uniformDistr (180 / 960) (200 / 960)))])
  ]

procTimesConst :: ConfigProcTimes
procTimesConst = ConfigProcTimes "Proc Times diverging 2 stages: Const" $
  map (second $ filter ((`elem` productTypes) . fst)) $ filter ((`elem` allBlocks) . fst)
    [ (Machine 1, [(Product 1, return . const (timeFromDouble (100 / 960))), (Product 2, return . const (timeFromDouble (100 / 960)))])
    , (Machine 2, [(Product 1, return . const (timeFromDouble (150 / 960)))])
    , (Machine 3, [(Product 2, return . const (timeFromDouble (190 / 960)))])
    ]

ext: hs

max_stars_repo_name: schnecki/borl-releaser

max_issues_repo_name: schnecki/borl-releaser

max_forks_repo_name: schnecki/borl-releaser


text:
-- |
-- Module      : Jeopardy.Controller
-- Description : MVC for our game
-- Copyright   : (c) Jonatan H Sundqvist, year
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created date year

-- TODO | -
--        -

-- SPEC | -
--        -

module Jeopardy.Controller where

---------------------------------------------------------------------------------------------------
-- We'll need these
---------------------------------------------------------------------------------------------------
import Graphics.UI.Gtk                   --
import Graphics.Rendering.Cairo (liftIO, fill) --

import Data.IORef   --
import Data.Complex --
import Data.Maybe (listToMaybe, maybe)
import Data.List (findIndex)

import Text.Printf

import qualified Southpaw.Picasso.Palette as Palette
import qualified Southpaw.Interactive.Application as App

import Jeopardy.Graphics -- TODO: Better name
import Jeopardy.Core     -- TODO: Better name
import Jeopardy.Curator  --
import Tiler             --

---------------------------------------------------------------------------------------------------
-- Types
---------------------------------------------------------------------------------------------------
-- |
-- TODO: Input state (?)
-- data App      = App { _window :: Window, _canvas :: DrawingArea, _size :: (Int, Int), _state :: IORef AppState } --
data AppState = AppState { _game :: Game, _selected :: Maybe Int, _path :: [Complex Double] }                    --

---------------------------------------------------------------------------------------------------
-- Data
---------------------------------------------------------------------------------------------------
fps = 30 :: Int

---------------------------------------------------------------------------------------------------
-- Functions
---------------------------------------------------------------------------------------------------
-- Events -----------------------------------------------------------------------------------------
-- |
-- ondelete :: IO Bool
ondelete = do
    -- TODO: Uhmmm... what?
    liftIO $ do
	    mainQuit
	    putStrLn "Goodbye!"
	    return False

-- |
onmousemotion stateref = do
	(mx, my) <- eventCoordinates
	liftIO $ do
		let cursor = mx:+my --
		let radius = 18     --
		appstate@(AppState { _path=path }) <- readIORef stateref
		writeIORef stateref $ appstate { _selected=findIndex (within radius cursor) path }
	return False
	where within r p = (< r) . realPart . abs . subtract p

-- |
ondraw appstate = do
	--
	renderGame $ _game appstate
	-- Testing an unrelated tiling function
	renderPathWithJoints Palette.chartreuse Palette.darkviolet 18 5 $ _path appstate
	assuming (_selected appstate) $ \ sel -> do
		renderCircle 22 $ _path appstate !! sel
		Palette.choose Palette.limegreen
		fill

		renderCircle 22 $ _path appstate !! opposite sel
		Palette.choose Palette.limegreen
		fill
	where assuming (Just a) f = f a
	      assuming  _       _ = return ()
	      opposite n          = (n + 3) `mod` (length $ _path appstate)

-- |
onanimate :: DrawingArea -> IORef AppState -> IO Bool
onanimate canvas stateref = do
    widgetQueueDraw canvas
    return True

---------------------------------------------------------------------------------------------------
-- |
createApp :: IO (App.App)
createApp = App.createWindowWithCanvas 650 650 $ AppState { _game=createGame,
		                                                    _selected=Nothing,
		                                                    _path=[ (200:+200)+(70:+0)*(cos θ:+sin θ) | θ <- [0, τ/6..(τ*5/6)] ] }

---------------------------------------------------------------------------------------------------
-- |
-- TODO: Rename, move (?)
mainGTK :: IO ()
mainGTK = do
	(App.App { App._window=window, App._canvas=canvas, App._size=size, App._state=stateref }) <- createApp

	timeoutAdd (onanimate canvas stateref) (1000 `div` fps) >> return ()

	-- Events
	canvas `on` draw $ (liftIO $ readIORef stateref) >>= ondraw
	canvas `on` motionNotifyEvent $ onmousemotion stateref

	-- canvas `on` buttonPressEvent   $ onbuttonpress worldref
	-- canvas `on` buttonReleaseEvent $ onbuttonreleased worldref

	-- window `on` configureEvent $ onresize window worldref
	window `on` deleteEvent	$ ondelete
	-- window `on` keyPressEvent  $ onkeypress worldref

	mainGUI

max_stars_repo_name: SwiftsNamesake/Leopardy

max_issues_repo_name: SwiftsNamesake/Leopardy

max_forks_repo_name: SwiftsNamesake/Leopardy


text:
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE TypeFamilies #-}

-- |
-- Module      :  Spiral.SPL
-- Copyright   :  (c) 2016-2017 Drexel University
-- License     :  BSD-style
-- Maintainer  :  mainland@drexel.edu

module Spiral.SPL (
    module Spiral.Permutation,

    SPL(..),
    matrix,
    fromLists,
    fromFunction,

    diag,
    circ,
    skew,
    toep,
    permute,
    backpermute,

    extent,

    mXv,
    (#>),

    toMatrix,
    col,
    row,

    transpose,

    inverse,

    (<->),
    (<|>),
    block,

    (⊗),
    (×),
    (⊕),

    Z(..),
    (:.)(..),

    DIM0,
    DIM1,
    DIM2,
    ix1,
    ix2
  ) where

import Data.Complex
import Data.Monoid ((<>))
import Data.Typeable (Typeable)
import qualified Data.Vector as V
import Text.PrettyPrint.Mainland hiding ((<|>))
import Text.PrettyPrint.Mainland.Class

import qualified Spiral.Array as A
import Spiral.Array (D,
                     IArray,
                     M,
                     Matrix,
                     Vector,
                     manifest)
import qualified Spiral.Array.Operators.Mapping as A
import qualified Spiral.Array.Operators.Matrix as A
import Spiral.Array.Repr.Complex
import Spiral.Array.Shape
import Spiral.Exp
import Spiral.Permutation hiding (Inv)
import Spiral.RootOfUnity
import Spiral.Util.Pretty

-- | Wrap an array to ensure it is unconditionally Show-able.
newtype ShowArray r sh e = ShowArray { unShowArray :: Array r sh e }

instance (Show sh, Show e, IArray r sh e) => Show (ShowArray r sh e) where
    show arr = show (manifest (unShowArray arr))

-- | An SPL transform. The type index is the type of the scalar values in the
-- transformed vector.
data SPL a where
    -- | An "embedded" array with an unknown representation.
    E :: IArray r DIM2 e => ShowArray r DIM2 e -> SPL e

    -- | The $n \times n$ identity matrix.
    I :: Num e => Int -> SPL e

    -- | Transpose
    T :: SPL e -> SPL e

    -- | Inverse
    Inv :: (Eq e, Fractional e) => SPL e -> SPL e

    -- | A permutation
    Pi :: Permutation -> SPL e

    -- | The rotation matrix
    Rot :: Floating a => a -> SPL a

    -- | A diagonal matrix
    Diag :: V.Vector e -> SPL e

    -- | The $n \times n$ diagonal matrix with constant diagonal elements.
    KDiag :: Int -> e -> SPL e

    -- | Vertical stacking of transforms
    Above :: SPL e -> SPL e -> SPL e

    -- | Horizontal stacking of transforms
    Beside :: SPL e -> SPL e -> SPL e

    -- | Kronecker product
    Kron :: SPL e -> SPL e -> SPL e

    -- | Direct sum
    DSum :: SPL e -> SPL e -> SPL e

    -- | Matrix product
    Prod :: SPL e -> SPL e -> SPL e

    -- | Circulant matrix given first column
    Circ :: V.Vector e -> SPL e

    -- | Skew-Circulant matrix given first column
    Skew :: V.Vector e -> SPL e

    -- | Toeplitz matrix
    Toep :: V.Vector e -> SPL e

    -- | Make a complex transform into a real transform.
    Re :: RealFloatConst a => SPL (Exp (Complex a)) -> SPL (Exp a)

    -- | The 2x2 DFT
    F2 :: SPL e

    -- | The nxn DFT matrix
    DFT :: RootOfUnity a => Int -> SPL a

    -- | The nxn Inverse DFT matrix
    DFT' :: RootOfUnity a => Int -> SPL a

    -- | The nxn "rotated" DFT matrix
    F :: RootOfUnity a => Int -> a -> SPL a

    -- | The nxn "rotated" inverse DFT matrix
    F' :: RootOfUnity a => Int -> a -> SPL a

deriving instance Show e => Show (SPL e)
deriving instance Typeable e => Typeable (SPL e)

-- | Embed any 'Matrix' as an SPL term.
matrix :: IArray r DIM2 e
       => Array r DIM2 e
       -> SPL e
matrix = E . ShowArray

-- | Embed a matrix created from a list of lists of values.
fromLists :: [[e]] -> SPL e
fromLists = matrix . A.fromLists

-- | Embed a matrix created from a function mapping indices to elements.
fromFunction :: DIM2 -> (DIM2 -> e) -> SPL e
fromFunction sh = matrix . A.fromFunction sh

-- | Create a diagonal matrix
diag :: [a] -> SPL a
diag = Diag . V.fromList

-- | Create a circulant matrix from the first column
circ :: [a] -> SPL a
circ = Circ . V.fromList

-- | Create a skew circulant matrix from the first column
skew :: [a] -> SPL a
skew = Skew . V.fromList

-- | Create a Toeplitz matrix
toep :: [a] -> SPL a
toep = Toep . V.fromList

-- | Permute (scatter).
permute :: Permutation -> SPL e
permute = Pi

-- | Backpermute (gather).
backpermute :: Permutation -> SPL e
backpermute = Pi . invert

-- | Return the extent of an SPL transform.
extent :: SPL a -> DIM2
extent (E a)       = A.extent (unShowArray a)
extent (I n)       = ix2 n n
extent (T e)       = ix2 n m
                       where
                         Z :. m :. n = extent e
extent (Inv e)     = ix2 n n
                       where
                         Z :. n :. _ = extent e
extent (Pi p)      = ix2 (dim p) (dim p)
extent (Rot _)     = ix2 2 2
extent (Diag xs)   = ix2 n n
                       where
                         n = length xs
extent (KDiag n _) = ix2 n n

extent (Above a b) = ix2 (ra+rb) c
  where
    Z :. ra :. c  = extent a
    Z :. rb :. _c = extent b

extent (Beside a b) = ix2 r (ca+cb)
  where
    Z :. r  :. ca  = extent a
    Z :. _r :. cb = extent b

extent (Kron a b) = ix2 (m*p) (n*q)
  where
    Z :. m :. n = extent a
    Z :. p :. q = extent b

extent (DSum a b) = ix2 (m+p) (n+q)
  where
    Z :. m :. n = extent a
    Z :. p :. q = extent b

extent (Prod a b) = ix2 m q
  where
    Z :. m  :. _n = extent a
    Z :. _p :. q  = extent b

extent (Circ xs) = ix2 n n
  where
    n = length xs

extent (Skew xs) = ix2 n n
  where
    n = length xs

extent (Toep xs) = ix2 n n
  where
    n = (length xs + 1) `quot` 2

extent (Re a) = ix2 (2*m) (2*n)
  where
    Z :. m  :. n = extent a

extent F2 = ix2 2 2

extent (DFT n)  = ix2 n n
extent (DFT' n) = ix2 n n
extent (F n _)  = ix2 n n
extent (F' n _) = ix2 n n

-- | Transpose an SPL expression
transpose :: forall a . (Show a, Num a) => SPL a -> SPL a
-- transpose (E m)        = (matrix . A.transpose . unShowArray) m
transpose a@I{}        = a
transpose (T a)        = a
transpose (Inv a)      = Inv (transpose a)
transpose (Pi p)       = backpermute p
transpose (Rot alpha)  = Rot (-alpha)
transpose a@Diag{}     = a
transpose a@KDiag{}    = a
transpose (Above a b)  = Beside (transpose a) (transpose b)
transpose (Beside a b) = Above (transpose a) (transpose b)
transpose (Kron a b)   = Kron (transpose a) (transpose b)
transpose (DSum a b)   = DSum (transpose a) (transpose b)
transpose (Prod a b)   = Prod (transpose b) (transpose a)
transpose (Circ cs)    = circ (x:reverse xs)
                           where
                             x:xs = V.toList cs
transpose a@Skew{}     = a
transpose (Toep cs)    = Toep (V.reverse cs)
transpose (Re a)       = Re (transpose a)
transpose F2           = F2
transpose a@DFT{}      = a
transpose a@DFT'{}     = a
transpose a@F{}        = a
transpose a@F'{}       = a
transpose a            = T a

-- | Inverse of an SPL expression
inverse :: forall a . (Eq a, Fractional a) => SPL a -> SPL a
inverse a@I{}       = a
inverse (T a)       = T (inverse a)
inverse (Inv a)     = a
inverse (Pi p)      = backpermute p
inverse (Rot alpha) = Rot (-alpha)
inverse (Diag xs)   = Diag (fmap (1/) xs)
inverse (KDiag n k) = KDiag n (1/k)
inverse (Kron a b)  = Kron (inverse a) (inverse b)
inverse (Prod a b)  = Prod (inverse b) (inverse a)
inverse (Re a)      = Re (inverse a)
inverse F2          = F2
inverse (DFT n)     = DFT' n
inverse (DFT' n)    = DFT n
inverse (F n w)     = F' n w
inverse (F' n w)    = F n w
inverse a           = Inv a

-- | Compute the product, @A x@.
mXv :: forall r a .
       ( Num a
       , IArray r DIM1 a
       )
    => SPL a      -- ^ The SPL transform, @A@
    -> Vector r a -- ^ The vector, @x@
    -> Vector D a
mXv I{} x =
    A.delay x

mXv (Pi p) x =
    A.fromFunction (ix1 n) f
  where
    Z :. n = A.extent x

    f (Z :. i) = x A.! reindex i
      where
        reindex = toIdxMapping (invert p)

mXv (Diag d) x =
    A.fromVector d A..* x

mXv (KDiag n k) x =
    A.fromFunction (ix1 n) f
  where
    f (Z :. i) = k * x A.! i

mXv (Kron (I m) a) x =
    A.delay $
    A.concat [a #> A.slice (i*n') 1 n' x | i <- [0..m-1]]
  where
    Z :. _n :. n' = extent a

mXv (Kron a (I n)) x =
    permute (L (m*n) m) #> Kron (I n) a #> permute (L (m'*n) n) #> x
  where
    Z :. m :. m' = extent a

mXv (Kron a b) x =
    (I m ⊗ b) #> (a ⊗ I n') #> x
  where
    Z :. m :. _  = extent a
    Z :. _ :. n' = extent b

mXv (DSum a b) x =
    A.delay $
    (a #> A.slice 0 1 m x) A.++ (b #> A.slice m 1 n x)
  where
    Z :. _ :. m = extent a
    Z :. _ :. n = extent b

mXv (Prod a b) x =
    a #> b #> x

mXv a x =
    toMatrix a `A.mXv` x

-- | Compute the matrix-vector product, @A x@.
infixr 8 #>
(#>) :: forall r a .
        ( Num a
        , IArray r DIM1 a
        )
     => SPL a
     -> Vector r a
     -> Vector D a
(#>) = mXv

-- | Convert an SPL transform to an explicit matrix.
toMatrix :: forall e . Num e => SPL e -> Matrix M e
toMatrix (E a) =
    manifest (unShowArray a)

toMatrix (T e) = manifest $ A.transpose $ toMatrix e

toMatrix (Inv e) =
    case A.inverse $ toMatrix e of
      Nothing -> error "Matrix has no inverse"
      Just a  -> A.manifest a

toMatrix (Diag xs) =
    manifest $ A.fromFunction (ix2 n n) f
  where
    n = length xs

    f (Z :. i :. j) | i == j    = xs V.! i
                    | otherwise = 0

toMatrix (KDiag n e) =
    manifest $ A.fromFunction (ix2 n n) f
  where
    f (Z :. i :. j) | i == j    = e
                    | otherwise = 0

toMatrix (Above a b) =
    A.manifest $ toMatrix a A.<-> toMatrix b

toMatrix (Beside a b) =
    A.manifest $ toMatrix a A.<|> toMatrix b

toMatrix (Kron a b) =
    A.manifest $ A.kronecker (toMatrix a) (toMatrix b)

toMatrix (DSum a b) =
    A.manifest $ A.directSum (toMatrix a) (toMatrix b)

toMatrix (Prod a b) =
    A.manifest $ A.mXm (toMatrix a) (toMatrix b)

toMatrix (Circ xs) =
    manifest $ A.fromFunction (ix2 n n) f
  where
    n = length xs

    f (Z :. i :. j) = xs V.! ((i-j) `mod` n)

toMatrix (Skew xs) =
    manifest $ A.fromFunction (ix2 n n) f
  where
    n = length xs

    f (Z :. i :. j) = xs V.! ((i+j) `mod` n)

toMatrix (Toep xs) =
    manifest $ A.fromFunction (ix2 n n) f
  where
    n = (length xs + 1) `quot` 2

    f (Z :. i :. j) = xs V.! (i-j+n-1)

toMatrix (I n) =
    manifest $ A.fromFunction (ix2 n n) f
  where
    f (Z :. i :. j) | i == j    = 1
                    | otherwise = 0

toMatrix (Rot alpha) =
    A.matrix [[cos alpha, -(sin alpha)],
              [sin alpha, cos alpha]]

toMatrix (Pi p) =
    manifest $ A.fromFunction (ix2 (dim p) (dim p)) f
  where
    f (Z :. i :. j) | g j == i  = 1
                    | otherwise = 0

    g = toIdxMapping p

toMatrix (Re a) = manifest (RE (toMatrix a))

toMatrix F2 =
    A.matrix [[1,  1],
              [1, -1]]

toMatrix (DFT n)  = toMatrix (F n (omega n))
toMatrix (DFT' n) = toMatrix (F' n (omega n))

toMatrix (F n w) = manifest $ A.fromFunction (ix2 n n) f
  where
    f (Z :. i :. j) = w ^ (i*j)

toMatrix (F' n w) = toMatrix (KDiag n (1/fromIntegral n) × F n (1/w))

pprArgs :: Pretty a => [a] -> Doc
pprArgs = parens . commasep . map ppr

instance (Num e, Pretty e) => Pretty (SPL e) where
    pprPrec p (E a)        = pprPrec p (manifest (unShowArray a))
    pprPrec _ (I n)        = text "I_" <> ppr n
    pprPrec _ (T e)        = pprPrec 10 e <> text "^T"
    pprPrec _ (Inv e)      = pprPrec 10 e <> text "^-1"
    pprPrec _ (Pi p)       = ppr p
    pprPrec _ (Rot alpha)  = text "R_" <> ppr alpha
    pprPrec _ (Diag xs)    = text "diag" <> pprArgs (V.toList xs)
    pprPrec _ (KDiag n e)  = text "kdiag" <> parens (commasep [ppr n, ppr e])
    pprPrec p (Above a b)  = infixop p AboveOp a b
    pprPrec p (Beside a b) = infixop p BesideOp a b
    pprPrec p (Kron a b)   = infixop p KOp a b
    pprPrec p (DSum a b)   = infixop p DSOp a b
    pprPrec p (Prod a b)   = infixop p POp a b
    pprPrec _ (Circ xs)    = text "circ" <> pprArgs (V.toList xs)
    pprPrec _ (Skew xs)    = text "skew" <> pprArgs (V.toList xs)
    pprPrec _ (Toep xs)    = text "toep" <> pprArgs (V.toList xs)
    pprPrec _ (Re a)       = text "Re" <> parens (ppr a)
    pprPrec _ F2           = text "F_2"
    pprPrec _ (DFT n)      = text "DFT_" <> ppr n
    pprPrec _ (DFT' n)     = text "DFT'_" <> ppr n
    pprPrec _ (F n w)      = text "F_" <> ppr n <> parens (ppr w)
    pprPrec _ (F' n w)     = text "F'_" <> ppr n <> parens (ppr w)

data MatrixBinop = AboveOp
                 | BesideOp
                 | KOp
                 | DSOp
                 | POp
  deriving (Eq, Ord, Show)

instance HasFixity MatrixBinop where
    fixity AboveOp  = infixl_ 7
    fixity BesideOp = infixl_ 7
    fixity KOp      = infixl_ 7
    fixity DSOp     = infixl_ 6
    fixity POp      = infixl_ 7

instance Pretty MatrixBinop where
    ppr AboveOp  = text "<->"
    ppr BesideOp = text "<|>"
    ppr KOp      = char '⊗'
    ppr DSOp     = char '⊕'
    ppr POp      = char '×'

-- | Extract a row of a matrix
row :: Matrix M e -> Int -> V.Vector e
row e i = V.slice (i*n) n xs
  where
    M (Z :. _m :. n) xs = manifest e

-- | Extract a column of a matrix
col :: forall e . Matrix M e -> Int -> V.Vector e
col e j = V.generate m f
  where
    M (Z :. m :. n) xs = manifest e

    f :: Int -> e
    f i = xs V.! (i*n + j)

-- | Vertical stacking of transforms
infixr 7 <->
(<->) :: SPL e -> SPL e -> SPL e
(<->) = Above

-- | Horizontal stacking of transforms
infixr 7 <|>
(<|>) :: SPL e -> SPL e -> SPL e
(<|>) = Beside

-- | Create a block matrix
block :: SPL e -> SPL e -> SPL e -> SPL e -> SPL e
block tl tr bl br = (tl <|> tr) <-> (bl <|> br)

-- | Alias for Kronecker product.
infixl 7 ⊗
(⊗) :: SPL e -> SPL e -> SPL e
I m ⊗ (Kron (I n) y) = I (m*n) ⊗ y
I 1 ⊗ y              = y
x   ⊗ I 1            = x
x   ⊗ y              = Kron x y

-- | Alias for matrix direct sum.
infixl 6 ⊕
(⊕) :: SPL e -> SPL e -> SPL e
I m     ⊕ I n     = I (n + m)
Diag xs ⊕ Diag ys = Diag (xs <> ys)
x       ⊕ y       = DSum x y

-- | Alias for matrix product.
infixl 7 ×
(×) :: SPL e -> SPL e -> SPL e
(×) = Prod

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
{-# LANGUAGE AllowAmbiguousTypes, DataKinds, RankNTypes, TypeFamilies,
             TypeOperators #-}
{-# OPTIONS_GHC -Wno-missing-export-lists #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}
-- | Shaped tensor-based implementation of Recurrent Neural Network
-- for classification of MNIST digits. Sports 2 hidden layers.
module HordeAd.Tool.MnistRnnShaped where

import Prelude

import           Control.Monad (foldM)
import qualified Data.Array.DynamicS as OT
import           Data.Array.Internal (valueOf)
import qualified Data.Array.Shape
import qualified Data.Array.Shaped as OSB
import qualified Data.Array.ShapedS as OS
import           Data.Proxy (Proxy)
import qualified Data.Vector.Generic as V
import           GHC.TypeLits (KnownNat)
import           Numeric.LinearAlgebra (Vector)
import qualified Numeric.LinearAlgebra as HM

-- until stylish-haskell accepts NoStarIsType
import qualified GHC.TypeLits

import HordeAd.Core.DualNumber
import HordeAd.Core.Engine
import HordeAd.Core.PairOfVectors (DualNumberVariables, varS)
import HordeAd.Tool.MnistData

zeroStateS
  :: (IsScalar d r, OS.Shape sh)
  => (DualNumber d (OS.Array sh r)  -- state
      -> a)
  -> a
zeroStateS f = f (konstS 0)

unrollLastS :: forall state c w m d r k rest.
               (DualMonad d r m, KnownNat k, OS.Shape rest)
            => (state -> OS.Array rest r -> w -> m (c, state))
            -> (state -> OS.Array (k : rest) r -> w -> m (c, state))
unrollLastS f s0 xs w = do
  let g :: (c, state) -> OS.Array rest r -> m (c, state)
      g (_, s) x = f s x w
  foldM g (undefined, s0) $ OSB.toList $ OS.unravel xs

type LayerWeigthsRNN in_width out_width d r =
  ( DualNumber d (OS.Array '[out_width, in_width] r)   -- input weight
  , DualNumber d (OS.Array '[out_width, out_width] r)  -- state weight
  , DualNumber d (OS.Array '[out_width] r) )           -- bias

rnnMnistLayerS
  :: forall in_width out_width batch_size d r m.
     ( DualMonad d r m
     , KnownNat in_width, KnownNat out_width, KnownNat batch_size )
  => DualNumber d (OS.Array '[out_width, batch_size] r)  -- in state
  -> DualNumber d (OS.Array '[in_width, batch_size] r)  -- in
  -> LayerWeigthsRNN in_width out_width d r
  -> m ( DualNumber d (OS.Array '[out_width, batch_size] r)  -- out
       , DualNumber d (OS.Array '[out_width, batch_size] r) )  -- out state
rnnMnistLayerS s x (wX, wS, b) = do
  let y = wX <>$ x + wS <>$ s + asColumnS b
  yTanh <- returnLet $ tanh y
  return (yTanh, yTanh)

rnnMnistTwoS
  :: forall out_width batch_size d r m.
     (DualMonad d r m, KnownNat out_width, KnownNat batch_size)
  => DualNumber d (OS.Array '[2 GHC.TypeLits.* out_width, batch_size] r)
       -- initial state
  -> OS.Array '[SizeMnistWidth, batch_size] r
  -> ( LayerWeigthsRNN SizeMnistWidth out_width d r
     , LayerWeigthsRNN out_width out_width d r )
  -> m ( DualNumber d (OS.Array '[out_width, batch_size] r)
       , DualNumber d (OS.Array '[2 GHC.TypeLits.* out_width, batch_size] r) )
           -- final state
rnnMnistTwoS s x ((wX, wS, b), (wX2, wS2, b2)) = do
  let s1 = sliceS @0 @out_width s
      s2 = sliceS @out_width @out_width s
  (vec1, s1') <- rnnMnistLayerS s1 (constant x) (wX, wS, b)
  (vec2, s2') <- rnnMnistLayerS s2 vec1 (wX2, wS2, b2)
  s3 <- returnLet $ appendS s1' s2'
  return (vec2, s3)

rnnMnistZeroS
  :: forall out_width batch_size d r m.
     (DualMonad d r m, KnownNat out_width, KnownNat batch_size)
  => OS.Array '[SizeMnistHeight, SizeMnistWidth, batch_size] r
  -- All below is the type of all paramters of this nn. The same is reflected
  -- in the length function below and read from variables further down.
  -> ( LayerWeigthsRNN SizeMnistWidth out_width d r
     , LayerWeigthsRNN out_width out_width d r )
  -> DualNumber d (OS.Array '[SizeMnistLabel, out_width] r)
  -> DualNumber d (OS.Array '[SizeMnistLabel] r)
  -> m (DualNumber d (OS.Array '[SizeMnistLabel, batch_size] r))
rnnMnistZeroS xs ((wX, wS, b), (wX2, wS2, b2)) w3 b3 = do
  (out, _s) <- zeroStateS (unrollLastS rnnMnistTwoS) xs
                          ((wX, wS, b), (wX2, wS2, b2))
  returnLet $ w3 <>$ out + asColumnS b3

rnnMnistLenS
  :: forall out_width. KnownNat out_width
  => Proxy out_width -> (Int, [Int], [(Int, Int)], [OT.ShapeL])
rnnMnistLenS _ =
  ( 0
  , []
  , []
  , [ Data.Array.Shape.shapeT @'[out_width, SizeMnistWidth]
    , Data.Array.Shape.shapeT @'[out_width, out_width]
    , Data.Array.Shape.shapeT @'[out_width]
    , Data.Array.Shape.shapeT @'[out_width, out_width]
    , Data.Array.Shape.shapeT @'[out_width, out_width]
    , Data.Array.Shape.shapeT @'[out_width]
    , Data.Array.Shape.shapeT @'[SizeMnistLabel, out_width]
    , Data.Array.Shape.shapeT @'[SizeMnistLabel]
    ]
  )

rnnMnistS
  :: forall out_width batch_size d r m.
     (DualMonad d r m, KnownNat out_width, KnownNat batch_size)
  => OS.Array '[SizeMnistHeight, SizeMnistWidth, batch_size] r
  -> DualNumberVariables d r
  -> m (DualNumber d (OS.Array '[SizeMnistLabel, batch_size] r))
rnnMnistS xs variables = do
  let wX = varS variables 0
      wS = varS variables 1
      b = varS variables 2
      wX2 = varS variables 3
      wS2 = varS variables 4
      b2 = varS variables 5
      w3 = varS variables 6
      b3 = varS variables 7
  rnnMnistZeroS @out_width xs ((wX, wS, b), (wX2, wS2, b2)) w3 b3

rnnMnistLossFusedS
  :: forall out_width batch_size d r m.
     (DualMonad d r m, KnownNat out_width, KnownNat batch_size)
  => Proxy out_width
  -> MnistDataBatchS batch_size r
  -> DualNumberVariables d r
  -> m (DualNumber d r)
rnnMnistLossFusedS _ (glyphS, labelS) variables = do
  let xs = OS.transpose @'[2, 1, 0] glyphS
  result <- rnnMnistS @out_width xs variables
  let targets2 = HM.tr $ HM.reshape (valueOf @SizeMnistLabel)
                       $ OS.toVector labelS
  vec <- lossSoftMaxCrossEntropyL targets2 (fromS2 result)
  returnLet $ scale (recip $ fromIntegral (valueOf @batch_size :: Int))
            $ sumElements0 vec

rnnMnistTestS
  :: forall out_width batch_size r.
     (IsScalar 'DModeGradient r, KnownNat out_width, KnownNat batch_size)
  => Proxy r -> Proxy out_width
  -> MnistDataBatchS batch_size r
  -> Domains r
  -> r
rnnMnistTestS _ _ (glyphS, labelS) parameters =
  let xs = OS.transpose @'[2, 1, 0] glyphS
      outputS = primalValue (rnnMnistS @out_width xs) parameters
      outputs = map OS.toVector $ OSB.toList $ OS.unravel
                $ OS.transpose @'[1, 0] $ outputS
      labels = map OS.toVector $ OSB.toList $ OS.unravel labelS
      matchesLabels :: Vector r -> Vector r -> Int
      matchesLabels output label | V.maxIndex output == V.maxIndex label = 1
                                 | otherwise = 0
  in fromIntegral (sum (zipWith matchesLabels outputs labels))
     / fromIntegral (valueOf @batch_size :: Int)

max_stars_repo_name: Mikolaj/horde-ad

max_issues_repo_name: Mikolaj/horde-ad

max_forks_repo_name: Mikolaj/horde-ad


text:
{-# LANGUAGE CPP #-}

-- | Binary lists are lists whose number of elements is a power of two.
--   This data structure is efficient for some computations like:
--
-- * Splitting a list in half.
--
-- * Appending two lists of the same length.
--
-- * Extracting an element from the list.
--
--   All the functions exported are total except for 'fromListWithDefault'.
--   It is impossible for the user of this library to create a binary list
--   whose length is /not/ a power of two.
--
--   Since many names in this module clash with the names of some "Prelude"
--   functions, you probably want to import this module this way:
--
-- > import Data.BinaryList (BinList,Exponent)
-- > import qualified Data.BinaryList as BL
--
--   Remember that binary lists are an instance of the 'Foldable' and 'Traversable'
--   classes. If you are missing a function here, look for functions using those
--   instances.
--
--   Note that some functions like 'replicate', 'generate', or 'take', don't use
--   the length of the list as argument, but the exponent of its length expressed
--   as a power of two. Throughout this document, this is referred
--   as the /length exponent/. For example, if the list has length 16, its length exponent
--   is 4 since 2^4 = 16. Therefore @replicate 4 0@ will create a list with 16 zeroes.
--   Keep this in mind when using this library. Note as well that this implies that
--   there is no need to check that the length argument is or is not a power of two.
--
module Data.BinaryList (
    -- * Type
    BinList
  , Exponent
    -- * Construction
  , singleton
  , append
  , replicate
  , replicateA
  , replicateAR
  , generate
  , generateM
    -- * Queries
  , lengthExponent
#if !MIN_VERSION_base(4,8,0)
  , length
#endif
  , lookup
  , head
  , last
    -- * Deconstruction
  , split
  , take
  , takeEnd
    -- * Transformation
  , replace
  , reverse
    -- * Tuples
  , joinPairs
  , disjoinPairs
    -- * Zipping and Unzipping
  , zip , unzip
  , zipWith
    -- * Lists
    -- ** From list
  , fromList
  , fromListWithDefault
  , fromListSplit
    -- ** To list
  , toListFilter
  , toListSegment
    -- * Others
  , traverseSegment
    -- * Example: Radix-2 FFT
    -- $fft
  ) where

import Prelude hiding
  ( map, head, last
  , reverse, replicate
  , take, lookup
  , zip, unzip, zipWith
  , foldr1, length, foldr
    )
#if MIN_VERSION_base(4,8,0)
import Prelude (length, foldr, foldr1)
#endif
import Data.Foldable (fold,toList)
import Foreign.Storable (sizeOf)
import Data.List (find)
import Data.BinaryList.Internal
import Control.Applicative.Backwards
import Control.Arrow ((***))
import Control.Monad.Trans.State
  ( StateT (..)
  , evalStateT ,evalState
  , runState
  , get ,modify )
import Control.Monad.Trans.Class (lift)
import Data.Functor.Identity (Identity (..))
import Control.Applicative.PhantomState
import Data.Word (Word8)

-- GHC-7.8 compatibility
#if !MIN_VERSION_base(4,8,0)

import qualified Prelude
import Data.Foldable (Foldable (..))
import Data.Traversable (Traversable (..))
import Control.Applicative (Applicative (..), (<$>))
import Data.Monoid (mappend)

-- | /O(1)/. Number of elements in the list.
length :: BinList a -> Int
length = (2^) . lengthExponent

#endif

-- | An exponent.
type Exponent = Word8

-- | /O(1)/. Build a list with a single element.
singleton :: a -> BinList a
singleton = ListEnd

-- | /O(1)/. Given a binary list @l@ with length @2^k@:
--
-- > lengthExponent l = k
--
lengthExponent :: BinList a -> Exponent
lengthExponent (ListNode n _ _) = n
lengthExponent (ListEnd _) = 0

{-# RULES
       "Data.BinaryList: length equality"
         forall xs ys . length xs == length ys = lengthExponent xs == lengthExponent ys
  #-}

-- | /O(log n)/. Lookup an element in the list by its index (starting from 0).
--   If the index is out of range, 'Nothing' is returned.
lookup :: BinList a -> Int -> Maybe a
lookup (ListNode n l r) i =
   let m = 2^(n-1) -- Number of elements in a single branch
   in  if i < m
          then lookup l i       -- Lookup in the left branch
          else lookup r $ i - m -- Lookup in the right branch
lookup (ListEnd x) 0 = Just x
lookup _ _ = Nothing

-- | /O(log n)/. Replace a single element in the list. If the index is
--   out of range, returns the original list.
replace :: Int -- ^ Index to look for
        -> a   -- ^ Element to insert
        -> BinList a -> BinList a
replace i0 y = go i0
  where
    go i (ListNode n l r) =
      let m = 2^(n-1)
      in  if i < m
             then ListNode n (go i l)          r
             else ListNode n       l (go (i-m) r)
    go 0 (ListEnd _) = ListEnd y
    go _ e = e

-- | /O(1)/. Append two binary lists. This is only possible
--   if both lists have the same length. If this condition
--   is not hold, 'Nothing' is returned.
append :: BinList a -> BinList a -> Maybe (BinList a)
append xs ys =
  let i = lengthExponent xs
  in  if i == lengthExponent ys
         then Just $ ListNode (i+1) xs ys
         else Nothing

-- | /O(1)/. Split a binary list into two sublists of half the length,
--   unless the list only contains one element. In that case, it
--   just returns that element.
split :: BinList a -> Either a (BinList a,BinList a)
split (ListNode _ l r) = Right (l,r)
split (ListEnd x) = Left x

-- | /O(log n)/. Calling @take n xs@ returns the first @min (2^n) (length xs)@ elements of @xs@.
take :: Exponent -> BinList a -> BinList a
take k xs@(ListNode n l _) = if k >= n then xs else take k l
take _ xs = xs

-- | /O(log n)/. Calling @takeEnd n xs@ returns the last @min (2^n) (length xs)@ elements of @xs@.
takeEnd :: Exponent -> BinList a -> BinList a
takeEnd k xs@(ListNode n _ r) = if k >= n then xs else takeEnd k r
takeEnd _ xs = xs

-- | Calling @replicateA n f@ builds a binary list collecting the results of
--   executing @2^n@ times the applicative action @f@.
replicateA :: Applicative f => Exponent -> f a -> f (BinList a)
replicateA n f = go n
  where
    go 0 = ListEnd <$> f
    go i = let b = go (i-1)
           in  ListNode <$> pure i <*> b <*> b

-- | The same as 'replicateA', but the actions are executed in reversed order.
replicateAR :: Applicative f => Exponent -> f a -> f (BinList a)
replicateAR n = forwards . replicateA n . Backwards

{-# RULES
      "Data.BinaryList: map reverse/replicateA"
         forall i f . map reverse (replicateA  i f) = replicateAR i f
  #-}

{-# RULES
      "Data.BinaryList: map reverse/replicateAR"
         forall i f . map reverse (replicateAR i f) = replicateA  i f
  #-}

-- | /O(log n)/. Calling @replicate n x@ builds a binary list with
--   @2^n@ occurences of @x@.
replicate :: Exponent -> a -> BinList a
replicate n = runIdentity . replicateA n . Identity

{-# RULES
      "Data.BinaryList: map/replicate"
         forall f n x . map f (replicate n x) = replicate n (f x)
  #-}

-- | /O(n)/. Build a binary list with the given length exponent (see 'lengthExponent')
--   by applying a function to each index.
generate :: Exponent -> (Int -> a) -> BinList a
generate l f = evalState (replicateA l $ fmap f get <* modify (+1)) 0

-- | Like 'generate', but the generator function returns a value in a 'Monad'.
--   Therefore, the result is as well contained in a 'Monad'.
generateM :: (Applicative m, Monad m) => Exponent -> (Int -> m a) -> m (BinList a)
generateM l f = evalStateT (replicateA l $ (get >>= lift . f) <* modify (+1)) 0

-- | /O(log n)/. Get the first element of a binary list.
head :: BinList a -> a
head (ListNode _ l _) = head l
head (ListEnd x) = x

-- | /O(log n)/. Get the last element of a binary list.
last :: BinList a -> a
last (ListNode _ _ r) = last r
last (ListEnd x) = x

-- | /O(n)/. Reverse a binary list.
reverse :: BinList a -> BinList a
{-# INLINE[2] reverse #-}
reverse (ListNode n l r) = ListNode n (reverse r) (reverse l)
reverse xs = xs

{-# RULES
      "Data.BinaryList: reverse/reverse"
         forall xs. reverse (reverse xs) = xs
  #-}

------------------------------
-- Transformations with tuples

-- | /O(n)/. Transform a list of pairs into a flat list. The
--   resulting list will have twice more elements than the
--   original.
joinPairs :: BinList (a,a) -> BinList a
{-# INLINE[1] joinPairs #-}
joinPairs (ListEnd (x,y)) = ListNode 1 (ListEnd x) (ListEnd y)
joinPairs (ListNode n l r) = ListNode (n+1) (joinPairs l) (joinPairs r)

-- | /O(n)/. Opposite transformation of 'joinPairs'. It halves
--   the number of elements of the input. As a result, when
--   applied to a binary list with a single element, it returns
--   'Nothing'.
disjoinPairs :: BinList a -> Maybe (BinList (a,a))
{-# INLINE [1] disjoinPairs #-}
disjoinPairs (ListEnd _) = Nothing
disjoinPairs xs = Just $ disjoinPairsNodes xs

disjoinPairsNodes :: BinList a -> BinList (a,a)
disjoinPairsNodes (ListNode _ (ListEnd x) (ListEnd y)) = ListEnd (x,y)
disjoinPairsNodes (ListNode n l r) = ListNode (n-1) (disjoinPairsNodes l) (disjoinPairsNodes r)
disjoinPairsNodes _ = error "disjoinPairsNodes: bug. Please, report this with an example input."

{-# RULES
      "Data.BinaryList: disjoinPairs/joinPairs"
         forall xs . disjoinPairs (joinPairs xs) = Just xs
  #-}

{-# RULES
      "Data.BinaryList: disjoinPairs/map/joinPairs"
         forall f xs . disjoinPairs (map f (joinPairs xs)) = Just (map (f *** f) xs)
  #-}

-- | Expression @pairBuilder f xs@ is equivalent to @joinPairs (map f xs)@, but does
--   not build any intermediate structure. Used for rewriting rules.
pairBuilder :: (a -> (b,b)) -> BinList a -> BinList b
{-# INLINE[0] pairBuilder #-}
pairBuilder f = go
  where
    go (ListEnd x) = let (a,b) = f x in ListNode 1 (ListEnd a) (ListEnd b)
    go (ListNode n l r) = ListNode (n+1) (go l) (go r)

{-# RULES
      "Data.BinaryList: joinPairs/map"
         forall f xs . joinPairs (map f xs) = pairBuilder f xs
  #-}

-- | Expression @zipAndJoing f g xs ys@ is equivalent to @pairBuilder f (zipWith g xs ys)@,
--   but does not build any intermediate structure. Used for rewriting rules.
zipAndJoin :: (c -> (d,d)) -> (a -> b -> c) -> BinList a -> BinList b -> BinList d
zipAndJoin f g = go
  where
    -- Recursion
    go xs@(ListNode n l r) ys@(ListNode n' l' r')
         -- If both lists have the same length, recurse assuming it
         -- to avoid comparisons.
       | n == n'   = ListNode (n+1) (goEquals l l') (goEquals r r')
         -- If the first list is larger, the second fits entirely in
         -- the left branch of the first.
       | n >  n'   = go l ys
         -- If the second list is larger, the first fits entirely in
         -- the left branch of the second.
       | otherwise = go xs l'
    go xs ys       = let (x,y) = f $ g (head xs) (head ys)
                     in  ListNode 1 (ListEnd x) (ListEnd y)
    -- Recursion assuming both lists have the same length
    goEquals (ListNode n l r) (ListNode _ l' r') =
                     ListNode (n+1) (goEquals l l') (goEquals r r')
    goEquals (ListEnd x) (ListEnd y) =
      let (x',y') = f $ g x y
      in  ListNode 1 (ListEnd x') (ListEnd y')
    goEquals _ _ = undefined -- This can't happen!

{-# RULES
      "Data.BinaryList: pairBuilder/zipWith"
         forall f g xs ys . pairBuilder f (zipWith g xs ys) = zipAndJoin f g xs ys
  #-}

------------------------
-- Zipping and Unzipping

-- | /O(n)/. Zip two binary lists using an operator.
zipWith :: (a -> b -> c) -> BinList a -> BinList b -> BinList c
{-# INLINE[1] zipWith #-}
zipWith f = go
  where
    -- Recursion
    go xs@(ListNode n l r) ys@(ListNode n' l' r')
         -- If both lists have the same length, recurse assuming it
         -- to avoid comparisons.
       | n == n'   = ListNode n (goEquals l l') (goEquals r r')
         -- If the first list is larger, the second fits entirely in
         -- the left branch of the first.
       | n >  n'   = go l ys
         -- If the second list is larger, the first fits entirely in
         -- the left branch of the second.
       | otherwise = go xs l'
    go xs ys       = ListEnd $ f (head xs) (head ys)
    -- Recursion assuming both lists have the same length
    goEquals (ListNode n l r) (ListNode _ l' r') =
                     ListNode n (goEquals l l') (goEquals r r')
    goEquals (ListEnd x) (ListEnd y) = ListEnd (f x y)
    goEquals _ _ = undefined -- this can't happen

-- | /O(n)/. Zip two binary lists in pairs.
zip :: BinList a -> BinList b -> BinList (a,b)
{-# INLINE zip #-}
zip = zipWith (,)

-- | /O(n)/. Unzip a binary list of pairs.
unzip :: BinList (a,b) -> (BinList a, BinList b)
{-# INLINE[1] unzip #-}
unzip (ListEnd (x,y)) = (ListEnd x, ListEnd y)
unzip (ListNode n l r) =
  let (la,lb) = unzip l
      (ra,rb) = unzip r
  in  (ListNode n la ra, ListNode n lb rb)

-- | Expression @unzipMap f xs@ is equivalent to @unzip (map f xs)@, but
--   does not create any intermediate structure.
unzipMap :: ((a,b) -> (c,d)) -> BinList (a,b) -> (BinList c,BinList d)
unzipMap f = go
  where
    go (ListEnd p) = ListEnd *** ListEnd $ f p
    go (ListNode n l r) =
      let (lc,ld) = go l
          (rc,rd) = go r
      in  (ListNode n lc rc, ListNode n ld rd)

{-# RULES
      "Data.BinaryList: unzip/map"
         forall f xs . unzip (map f xs) = unzipMap f xs
  #-}

-----------------------------
-- Transforming from/to lists

-- | /O(log n)/. Calculate the exponent of a positive integer number expressed
--   as a power of two.
exponentInBasisTwo :: Int -> Maybe Exponent
exponentInBasisTwo 1 = Just 0
exponentInBasisTwo n =
  if even n
     then fmap (+1) $ exponentInBasisTwo $ div n 2
     else Nothing

-- | /O(n)/. Build a binary list from a linked list. If the input list
--   has length different from a power of two, it returns 'Nothing'.
fromList :: [a] -> Maybe (BinList a)
fromList xs = fmap builder . exponentInBasisTwo $ Prelude.length xs
  where
    builder l = evalState (replicateA l $ StateT $ \(h:t) -> pure (h,t)) xs

-- | /O(1)/. This is the last exponent that has power of two defined in the type 'Int'.
--
-- /Note: This value is system dependent, since the type 'Int' varies in size/
-- /from system to system./
--
lastExponentOfTwo :: Exponent
lastExponentOfTwo = fromIntegral $ 8 * sizeOf (undefined :: Int) - 2

-- | /O(1)/. Calculate the next power of two exponent, if there is any. It is possible
--   to not find a next one since the type 'Int' is finite. If the input is
--   already a power of two, its exponent is returned.
nextExponentOfTwo :: Int -> Maybe Exponent
nextExponentOfTwo n = find (\i -> n <= 2^i) [0 .. lastExponentOfTwo]

-- | /O(n)/. Build a binary list from a linked list. If the input list
--   has length different from a power of two, fill to the next
--   power of two with a default element.
--
-- /Warning: this function crashes if the input list length is larger than any/
-- /power of two in the type 'Int'. However, this is very unlikely./
fromListWithDefault :: a -> [a] -> BinList a
fromListWithDefault e xs =
  let l = Prelude.length xs
  in  case nextExponentOfTwo l of
        Just n ->
          evalState (replicateA n $ StateT $
             \ys -> pure $ case ys of
                      (h:t) -> (h,t )
                      _     -> (e,[])
               ) xs
        _ -> error "[binary-list] fromListWithDefault: input list is too big."

-- | /O(n)/. Build a binary list from a linked list. It returns a binary list
--   with length @2 ^ n@ (where @n@ is the supplied 'Int' argument), and
--   the list of elements of the original list that were not used. If the
--   input list is shorter than @2 ^ n@, a default element will be used to
--   complete the binary list. This method for building binary lists is faster
--   than both 'fromList' and 'fromListWithDefault'.
fromListSplit :: a   -- ^ Default element
              -> Exponent -- ^ Length exponent
              -> [a] -- ^ Input list
              -> (BinList a, [a])
fromListSplit e n =
  runState $ replicateA n $ StateT $
    \xs -> pure $ case xs of
                    (h:t) -> (h,t )
                    _     -> (e,[])

-- | /O(n)/. Create a list from the elements of a binary list matching a given
--   condition.
toListFilter :: (a -> Bool) -> BinList a -> [a]
{-# INLINE toListFilter #-}
toListFilter c = foldr (\x -> if c x then (x:) else id) []

-- | /O(n)/. Create a list extracting a sublist of elements from a binary list.
toListSegment :: Int -> Int -> BinList a -> [a]
{-# INLINE toListSegment #-}
toListSegment s e xs = runPhantomState (traverseSegment (changeState . (:)) s e xs) []

-- | Apply an applicative action to every element in a segment of a binary list, from left to right.
traverseSegment :: Applicative f => (a -> f ()) -> Int -> Int -> BinList a -> f ()
{-# INLINE traverseSegment #-}
traverseSegment f s e xs
  | s > e = pure ()
  | e < 0 = pure ()
  | s >= length xs = pure ()
  | otherwise = traverseSegmentFromTo f (max 0 s) e xs

traverseSegmentFromTo :: Applicative f => (a -> f ()) -> Int -> Int -> BinList a -> f ()
{-# INLINE traverseSegmentFromTo #-}
traverseSegmentFromTo f = go
  where
    go s e (ListNode n l r) =
      let k = 2^(n-1)
      in  if s >= k
             -- Sublist is contained in right portion
             then go (s - k) (e - k) r
             else if e < k
                     -- Sublist is contained in left portion
                     then go s e l
                     -- Sublist is divided in both portions
                     else traverseSegmentFrom f s l *> traverseSegmentTo f (e - k) r
    go _ _ (ListEnd x) = f x

traverseSegmentFrom :: Applicative f => (a -> f ()) -> Int -> BinList a -> f ()
{-# INLINE traverseSegmentFrom #-}
traverseSegmentFrom f = go
  where
    go s (ListNode n l r) =
      let k = 2^(n-1)
      in  if s >= k
             -- Sublist is contained in right portion
             then go (s - k) r
             -- Sublist is divided in both portions, but right
             -- portion is taken entirely
             else go s l *> traverseFull f r
    go _ (ListEnd x) = f x

traverseSegmentTo :: Applicative f => (a -> f ()) -> Int -> BinList a -> f ()
{-# INLINE traverseSegmentTo #-}
traverseSegmentTo f = go
  where
    go e (ListNode n l r) =
      let k = 2^(n-1)
      in  if e < k
             -- Sublist is contained in left portion
             then go e l
             -- Sublist is divided in both portions, but left
             -- portion is taken entirely
             else traverseFull f l *> go (e - k) r
    go _ (ListEnd x) = f x

traverseFull :: Applicative f => (a -> f ()) -> BinList a -> f ()
{-# INLINE traverseFull #-}
traverseFull f = go
  where
    go (ListEnd x) = f x
    go (ListNode _ l r) = go l *> go r

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

-----------------------------
-- Show and Functor instances

instance Show a => Show (BinList a) where
  show = show . toList

{- Internal map

Although we encourage the use of 'fmap', we define fmap as a custom 'map'
function and inline 'fmap' to make them equivalent, so writing 'fmap' is
actually writing 'map'. We do this to use 'map' in rewriting rules.

-}

map :: (a -> b) -> BinList a -> BinList b
{-# INLINE[1] map #-}
map f = go
  where
    go (ListEnd x) = ListEnd (f x)
    go (ListNode n l r) = ListNode n (go l) (go r)

instance Functor BinList where
  {-# INLINE fmap #-}
  fmap = map

instance Foldable BinList where
  -- Folding
  foldr1 f = go
    where
      go (ListEnd x) = x
      go (ListNode _ l r) = f (go l) (go r)
  --
  fold = foldr1 mappend
  foldl1 = foldr1
  foldMap f = fold . fmap f
#if MIN_VERSION_base(4,8,0)
  length = (2^) . lengthExponent
#endif

instance Traversable BinList where
  sequenceA (ListEnd f) = ListEnd <$> f
  sequenceA (ListNode n l r) = ListNode <$> pure n <*> sequenceA l <*> sequenceA r

-----------------------------
-- Example: Radix-2 FFT

{- $fft

This is an example demonstrating the use of binary lists to calculate the Discrete
Fourier Transform of complex vectors with the Radix-2 Fast Fourier Transform algorithm.

> import Data.BinaryList (BinList)
> import qualified Data.BinaryList as BL
>
> import Data.Complex
> import Data.Maybe (fromJust)
>
> i :: Complex Double
> i = 0 :+ 1
>
> fft :: BinList (Complex Double) -> BinList (Complex Double)
> fft xs =
>   case BL.disjoinPairs xs of
>     Nothing -> xs
>     Just ps ->
>       let (evens,odds) = BL.unzip ps
>           n = BL.lengthExponent xs - 1
>           q = negate $ pi * i / fromIntegral (2^n)
>           twiddles = BL.generate n $ \k -> exp $ q * fromIntegral k
>           oddsfft = BL.zipWith (*) twiddles $ fft odds
>           evensfft = fft evens
>       in  fromJust $
>             BL.append (BL.zipWith (+) evensfft oddsfft)
>                       (BL.zipWith (-) evensfft oddsfft)

-}

max_stars_repo_name: Daniel-Diaz/binary-list

max_issues_repo_name: Daniel-Diaz/binary-list

max_forks_repo_name: Daniel-Diaz/binary-list


text:
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE GADTs                     #-}
{-# LANGUAGE MultiParamTypeClasses     #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE OverloadedStrings         #-}
{-# LANGUAGE PolyKinds                 #-}
{-# LANGUAGE RankNTypes                #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# OPTIONS_GHC -fwarn-incomplete-patterns #-}

module Math.Regression.LeastSquares where

import qualified Math.HMatrixUtils             as HU
import qualified Math.Regression.Regression    as RE

import qualified Polysemy                      as P
import qualified Knit.Effect.Logger            as Log

import qualified Data.List                     as List

import           Numeric.LinearAlgebra          ( (#>)
                                                , (<.>)
                                                , (<\>)
                                                )
import qualified Numeric.LinearAlgebra         as LA
import           Numeric.LinearAlgebra.Data     ( Matrix
                                                , R
                                                , Vector
                                                )

-- matrix dimensions given in (row x column) form
-- if A is (m x n), its transpose, A' is (n x m)
-- ||A|| is the Frobenius norm of A, the sum of the squares of the elements.
-- if A (m x l) and B (m x k) have the same number of rows, [A B] is the (m x (l+k)) matrix formed by appending the columns of B to the columns of A
-- A: an (m x n) matrix, of m observations of n observables
-- B: an (m x d) matrix, of observations (often d=1 and B is just a length m column-vector)
-- Q: an ((n+d)m x (n+d)m) matrix of covariances, the covariance matrix of [A B].  Given?  Computed?

-- we want an (n x d) matrix X which solves, or "best solves", AX = B.

-- OLS: minimize ||AX - B||
ordinaryLS
  :: Log.LogWithPrefixesLE effs
  => Bool
  -> Matrix R
  -> Vector R
  -> P.Sem effs (RE.RegressionResult R)
ordinaryLS withConstant mA vB = do
  let mAwc   = if withConstant then addBiasCol (LA.size vB) mA else mA -- add a constant, e.g., the b in y = mx + b
      (n, p) = LA.size mAwc
      dof    = n - p
  HU.checkVectorMatrix "b" "A" vB mAwc  -- vB <> mA is legal, b has same length as A has rows
  let vX  = mAwc <\> vB
      vU  = vB - (mAwc #> vX) -- residuals
      mse = (vU <.> vU) / realToFrac dof
      cov = LA.scale mse (LA.inv $ LA.tr mAwc LA.<> mAwc)
  RE.FitStatistics rSq aRSq fStat <- RE.goodnessOfFit p vB Nothing vU
  return $ RE.RegressionResult (RE.estimates cov vX)
                               (realToFrac dof)
                               mse
                               rSq
                               aRSq
                               fStat
                               cov

weightedLS
  :: Log.LogWithPrefixesLE effs
  => Bool
  -> Matrix R
  -> Vector R
  -> Vector R
  -> P.Sem effs (RE.RegressionResult R)
weightedLS withConstant mA vB vW = do
  HU.checkEqualVectors "b" "w" vB vW
  let mW     = LA.diag vW
      vWB    = mW #> vB
      mAwc   = if withConstant then addBiasCol (LA.size vB) mA else mA -- add a constant, e.g., the b in y = mx + b
      (_, p) = LA.size mAwc
  HU.checkVectorMatrix "b" "A" vB mAwc
  let mWA  = mW LA.<> mAwc
      vX   = mWA <\> vWB
      vU   = vB - (mAwc #> vX)
      vWU  = mW #> vU
      sumW = LA.sumElements vW
      effN = (sumW * sumW) / (vW <.> vW)
      mse  = effN / (effN - realToFrac p) * (vWU <.> vWU) / sumW
      cov  = LA.scale mse (LA.inv $ LA.tr mAwc LA.<> mAwc)
  RE.FitStatistics rSq aRSq fStat <- RE.goodnessOfFit p vB (Just vW) vU
  return $ RE.RegressionResult (RE.estimates cov vX)
                               (effN - realToFrac p)
                               mse
                               rSq
                               aRSq
                               fStat
                               cov

totalLS
  :: Log.LogWithPrefixesLE effs
  => Bool
  -> Matrix R
  -> Vector R
  -> P.Sem effs (RE.RegressionResult R)
totalLS withConstant mA vB = do
  let mAwc   = if withConstant then addBiasCol (LA.size vB) mA else mA -- add a constant, e.g., the b in y = mx + b
      (n, p) = LA.size mAwc
      dof    = realToFrac (n - p)
  HU.checkVectorMatrix "b" "Awc" vB mAwc
  let mAB      = mAwc LA.||| LA.asColumn vB
      (_, mV') = LA.rightSV mAB
--      gSV = sv V.! 0
--      tol = gSV * LA.eps
      sV22     = mV' LA.! p LA.! p
      vV12     = List.head $ LA.toColumns $ LA.subMatrix (0, p) (p, 1) mV' --LA.?? (LA.DropLast 1, LA.TakeLast 1)
      vX       = LA.scale (-1 / sV22) vV12 -- this is the TLS solution.  But in a shifted basis.??
--      mV2      = mV' LA.?? (LA.All, LA.TakeLast 1)
--      mABt     = mAB LA.<> mV2 LA.<> (LA.tr mV2)
--      mAt      = mABt LA.?? (LA.All, LA.DropLast 1)
      vBfit    = mAwc #> vX --(mAwc - mAt) #> vX -- this is confusing.
      vU       = vB - vBfit
--      (rSq, aRSq) = RE.goodnessOfFit p vB vU --(vB - mA #> vX)
      mse      = (vU <.> vU) / realToFrac (n - p)
      cov      = LA.scale mse (LA.inv $ LA.tr mAwc LA.<> mAwc)
  RE.FitStatistics rSq aRSq fStat <- RE.goodnessOfFit p vB Nothing vU
  return $ RE.RegressionResult (RE.estimates cov vX) dof mse rSq aRSq fStat cov

weightedTLS
  :: Log.LogWithPrefixesLE effs
  => Bool
  -> Matrix R
  -> Vector R
  -> Vector R
  -> P.Sem effs (RE.RegressionResult R)
weightedTLS withConstant mA vB vW = do
  HU.checkEqualVectors "b" "w" vB vW
  let mW     = LA.diag vW
      vWB    = mW #> vB
      mAwc   = if withConstant then addBiasCol (LA.size vB) mA else mA -- add a constant, e.g., the b in y = mx + b
      (_, p) = LA.size mAwc
  HU.checkVectorMatrix "b" "Awc" vB mAwc
  let mWA      = mW LA.<> mAwc
      mWAB     = mWA LA.||| LA.asColumn vWB
      (_, mV') = LA.rightSV mWAB
--      gSV = sv V.! 0
--      tol = gSV * LA.eps
      sV22     = mV' LA.! p LA.! p
      vV12     = List.head $ LA.toColumns $ LA.subMatrix (0, p) (p, 1) mV' --LA.?? (LA.DropLast 1, LA.TakeLast 1)
      vX       = LA.scale (-1 / sV22) vV12 -- this is the WTLS solution.  But in a shifted basis.??
--      mV2      = mV' LA.?? (LA.All, LA.TakeLast 1)
--      mWABt    = mWAB LA.<> mV2 LA.<> (LA.tr mV2)
--      mWAt     = mWABt LA.?? (LA.All, LA.DropLast 1)
--      vBfit    = mAwc #> vX -- (mWA - mWAt) #> vX -- this is confusing
      vU       = vB - (mAwc #> vX)
      vWU      = mW #> vU
      sumW     = LA.sumElements vW
      effN     = (sumW * sumW) / (vW <.> vW)
      dof      = effN - realToFrac p
      mse      = effN / dof * (vWU <.> vWU) / sumW
--      (rSq, aRSq) = RE.goodnessOfFit p vWB vWU --(vB - mA #> vX)
      cov      = LA.scale mse (LA.inv $ LA.tr mAwc LA.<> mAwc)
  RE.FitStatistics rSq aRSq fStat <- RE.goodnessOfFit p vB (Just vW) vU
  return $ RE.RegressionResult (RE.estimates cov vX) dof mse rSq aRSq fStat cov

addBiasCol :: Int -> Matrix R -> Matrix R
addBiasCol rows mA =
  let colList = List.replicate rows 1
  in  if LA.size mA == (0, 0)
        then LA.matrix 1 colList
        else LA.col colList LA.||| mA

max_stars_repo_name: teto/Frames-utils

max_issues_repo_name: teto/Frames-utils

max_forks_repo_name: teto/Frames-utils


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeFamilies #-}
module Numeric.Layer.Affine where

import Dependent.Size
import Numeric.Layer
import Numeric.Vector.Sized
import qualified Numeric.Vector.Sized as Sized
import GHC.TypeLits
import Data.Array.Accelerate.LLVM.Native

data Affine (n :: Nat) (m :: Nat) = Affine (SizedMatrix' n m) (SizedVector' m)
  deriving (Show)

instance (KnownNat n, KnownNat m) => Layer (Affine n m) where
  type Inputs (Affine n m) = ZZ '::. n
  type Outputs (Affine n m) = ZZ '::. m
  type Tape (Affine n m) = SizedArray (ZZ ::. n)
  type Gradient (Affine n m) = Affine n m
  forward (Affine ws b) x = (affine (use ws) (use b) x, x)
  backward (Affine ws b) x dy = (Affine dw db, dx)
    where
      dx = use ws #> dy
      dw = runSized run (dy `outer` x)
      db = runSized run (dy)
  applyGradient (Affine ws b) (Affine dw db)
    = Affine
      (runSized run $ Sized.zipWith (-) (use ws) (use dw))
      (runSized run $ Sized.zipWith (-) (use b) (use db))

affine
  :: (KnownNat n, KnownNat m)
  => SizedMatrix n m
  -> SizedVector m
  -> SizedVector n
  -> SizedVector m
affine ws b x = (x <# ws) .+. b

-- affine_dw dedy = dedw
-- affine_db dedy = dedb
-- affine_dx dedy = dedx

max_stars_repo_name: mixed-signals/mixed-signals

max_issues_repo_name: mixed-signals/mixed-signals

max_forks_repo_name: mixed-signals/mixed-signals


text:
{-# LANGUAGE FlexibleContexts,
             MultiParamTypeClasses #-}
module Network.Trainer.BackpropTrainer
( BackpropTrainer(..)

, backprop
, inputs
, outputs
, deltas
, hiddenDeltas
, calculateNablas
, fit
, evaluate
) where

import Network.Network
import Network.Network.FeedForwardNetwork
import Network.Neuron
import Network.Layer
import Network.Trainer

import Numeric.LinearAlgebra

-- | A BackpropTrainer performs simple backpropagation on a neural network.
--   It can be used as the basis for more complex trainers.
data BackpropTrainer = BackpropTrainer { eta :: Double
                                       , cost :: CostFunction
                                       , cost' :: CostFunction'
                                       }

-- | Declare the BackpropTrainer to be an instance of Trainer.
instance Trainer BackpropTrainer FeedForwardNetwork where
  fit s t n examples = foldl (backprop t) n $ s examples
  -- | Use the cost function to determine the error of a network
  evaluate t n example = (cost t) (snd example) (predict (fst example) n)

-- | Perform backpropagation on a single training data instance.
backprop :: BackpropTrainer -> FeedForwardNetwork -> [TrainingData] -> FeedForwardNetwork
backprop t n es =
  updateNetwork (length es) t (foldl (calculateNablas t n) emptyFeedForwardNetwork es) n

-- | Given the size of the minibatch, the trainer, the nablas for each layer, given
--   as a network, and the network itself, return a network with updated wieghts.
updateNetwork :: Int -> BackpropTrainer -> FeedForwardNetwork -> FeedForwardNetwork -> FeedForwardNetwork
updateNetwork mag t nablas n = addFeedForwardNetworks n
  (FeedForwardNetwork $ map (scaleLayer $ -1 * (eta t) / (fromIntegral mag)) (layers nablas))

-- | Calculate the nablas for a minibatch and return them as a network (so each
--   weight and bias gets its own nabla).
calculateNablas :: BackpropTrainer -> FeedForwardNetwork -> FeedForwardNetwork -> TrainingData -> FeedForwardNetwork
calculateNablas t n nablas e = FeedForwardNetwork $ map (updateLayer t) (zip3 (layers n) ds os)
  where ds = deltas t n e
        os = outputs (fst e) n

-- | The mapped function to update the weight and biases in a single layer
updateLayer :: BackpropTrainer -> (Layer, Vector Double, Vector Double) -> Layer
updateLayer t (l, delta, output) = Layer newWeight newBias n
  where n = neuron l
        newWeight = ((reshape 1 delta) <> (reshape (dim output) output))
        newBias = delta

-- | The outputs function scans over each layer of the network and stores the
--   activated results
outputs :: Vector Double -> FeedForwardNetwork -> [Vector Double]
outputs input network = scanl apply input (layers network)

-- | The inputs function performs a similar task to outputs, but returns a list
--   of vectors of unactivated inputs
inputs :: Vector Double -> FeedForwardNetwork -> [Vector Double]
inputs input network = if null (layers network) then []
  else unactivated : inputs activated (FeedForwardNetwork (tail $ layers network))
  where unactivated = weightMatrix layer <> input + biasVector layer
        layer = head $ layers network
        activated = mapVector (activation (neuron layer)) unactivated

-- | The deltas function returns a list of layer deltas.
deltas :: BackpropTrainer -> FeedForwardNetwork -> TrainingData -> [Vector Double]
deltas t n example = hiddenDeltas
  (FeedForwardNetwork (reverse (layers n))) outputDelta (tail $ reverse is)
    ++ [outputDelta]
  where outputDelta = costd (snd example) output *
          mapVector activationd lastInput
        costd = cost' t
        activationd = activation' (neuron (last (layers n)))
        output = last os
        lastInput = last is
        is = inputs (fst example) n
        os = outputs (fst example) n

-- | Compute the hidden layer deltas
hiddenDeltas :: FeedForwardNetwork -> Vector Double -> [Vector Double] -> [Vector Double]
hiddenDeltas n prevDelta is = if length (layers n) <= 1 then []
  else delta : hiddenDeltas rest delta (tail is)
  where rest = FeedForwardNetwork (tail $ layers n)
        delta = (trans w) <> prevDelta * spv
        w = weightMatrix (head $ layers n)
        spv = mapVector (activation' (neuron (head $ layers n))) (head is)

max_stars_repo_name: AkatsukiSirius/LambdaNet

max_issues_repo_name: world-admin/LambdaNet

max_forks_repo_name: world-admin/LambdaNet


text:
{-# LANGUAGE RecordWildCards #-}

{-|
Module      :  HABQTlib.Data.Particle

Data structures and functions that deal with storing and processing particle
hierarchies.

/Warning/: functions in this module assume that the 'ptsParticles' is non-empty
and 'NumberOfParticles', 'Dim', and 'Rank' are positive, no validation is
performed. If you use them directly, instead of employing API from
"HABQTlib", you must ensure those assumptions hold.
-}
module HABQTlib.Data.Particle
  ( Particles(..)
  , genParticles
  , updateParticles
  , ParticleHierarchy
  , initialiseParticleHierarchy
  , updateParticleHierarchy
  , getMixedEstimate
  , foldOverPts
  , reduceParticlesToMean
  , effectiveSize
  , ResampleArgs(..)
  , resampleMultinom
  , resample
  , ecdf
  , icdf
  , nudgeParticle
  ) where

import Control.Monad (when)
import Control.Newtype.Generics (over)
import Data.Bool.HT (if', select)
import qualified Data.Vector as V
import HABQTlib.Data
import HABQTlib.RandomStates
import Numeric.LinearAlgebra (Complex(..))
import qualified Numeric.LinearAlgebra as LA
import qualified System.Random.MWC as MWC
import Text.Printf (printf)

-- | A vector of weighed density matrices is stored along with their rank and
-- number. 'ptsWeight' corresponds to the collective weight of particles of
-- rank 'ptsRank' in the hierarchical model, it is not the sum of individual
-- weights of particles (that is normalised to unity after every update).
data Particles = Particles
  { ptsRank :: Rank
  , ptsWeight :: Weight
  , ptsNumber :: NumberOfParticles
  , ptsParticles :: V.Vector WeighedDensityMatrix
  } deriving (Show)

-- | Particle hierarchy is described by a vector of 'Particles'.
type ParticleHierarchy = V.Vector Particles

-- | Generates random particles (according to induced measure).
genParticles :: Dim -> Rank -> NumberOfParticles -> IO Particles
genParticles d r n =
  let w = 1 / fromIntegral n
   in Particles r 1 n . fmap (mkWDM w) <$> V.replicateM n (genDM d r)

-- | Summarise particles to a mean estimate (weighed by the corresponding
-- hierarchical weight of the rank).
reduceParticlesToMean :: Particles -> WeighedDensityMatrix
reduceParticlesToMean Particles {..} =
  let wdm = V.foldl1' (<+>) ptsParticles
      wdmw = over WeighedDensityMatrix (\(w, dm) -> (w * ptsWeight, dm)) wdm
   in truncateRank ptsRank wdmw

-- | Map density matrices, combine them with their weights, and then perform a
-- (strict left) fold of results.
foldOverPts ::
     (DensityMatrix -> a) -- ^ function to map over density matrices
  -> (Weight -> a -> b) -- ^ function to combine weights with results of mapping
  -> (c -> b -> c) -- ^ fold funciton
  -> c -- ^ seed value for folding
  -> Particles
  -> c
foldOverPts f wf fld z Particles {..} =
  let wm (WeighedDensityMatrix (w, dm)) = wf w (f dm)
   in V.foldl' (\l r -> fld l (wm r)) z ptsParticles

fullDataLogLikelihood :: [PureStateVector] -> DensityMatrix -> Double
fullDataLogLikelihood vs dm =
  let lps = map (log . (`pureStateLikelihood` dm)) vs
   in sum lps

-- | Given a measurement result, perform a Bayesian update over the particles.
updateParticles :: PureStateVector -> Particles -> Particles
updateParticles sv pts@Particles {..} =
  let updateF :: WeighedDensityMatrix -> WeighedDensityMatrix
      updateF (WeighedDensityMatrix (w, dm)) = WeighedDensityMatrix (wnew, dm)
        where
          wnew = w * pureStateLikelihood sv dm
      upts = V.map updateF ptsParticles
      uw = V.foldl' (\acc (WeighedDensityMatrix (w, _)) -> acc + w) 0 upts
      npts = V.map (over WeighedDensityMatrix (\(w, dm) -> (w / uw, dm))) upts
   in pts {ptsWeight = ptsWeight * uw, ptsParticles = npts}

-- | Helper function that generates random particles of each applicable rank.
initialiseParticleHierarchy :: Dim -> NumberOfParticles -> IO ParticleHierarchy
initialiseParticleHierarchy d n = V.generateM d (\r -> genParticles d (r + 1) n)

-- | Given a measurement result, update all particles, then normalise resulting
-- hierarchical weights to sum to unity.
updateParticleHierarchy ::
     PureStateVector -> ParticleHierarchy -> ParticleHierarchy
updateParticleHierarchy sv ph =
  let uph = V.map (updateParticles sv) ph
      wgts = V.map ptsWeight uph
      nwgts = V.map (/ V.sum wgts) wgts
   in V.zipWith (\x w -> x {ptsWeight = w}) uph nwgts

-- | Summarise the whole particle hierarchy into one mean Bayesian estimate.
getMixedEstimate :: ParticleHierarchy -> DensityMatrix
getMixedEstimate ph =
  let rankEstimates = V.map reduceParticlesToMean ph
      WeighedDensityMatrix (_, result) = V.foldl1' (<+>) rankEstimates
   in result

-- | Calculate the effective sample size of particles (weights don’t
-- necessarily have to be normalised).
effectiveSize :: Particles -> Double
effectiveSize Particles {..} =
  let ss = V.sum . V.map ((^ (2 :: Int)) . fst . getWDM) $ ptsParticles
      wa = V.foldl' (flip ((+) . fst . getWDM)) 0 ptsParticles
   in wa ^ (2 :: Int) / ss

-- | Nudges a particle by mixing the state together with some randomly
-- generated pure state. Relative weight of the random component determines how
-- “close” a nudged particle is to the original one.
nudgeParticle ::
     Dim
  -> Weight -- ^ Relative weight (from 0 to 1) of random component
  -> WeighedDensityMatrix
  -> IO WeighedDensityMatrix
nudgeParticle dim weightFraction (WeighedDensityMatrix (w, dm)) = do
  DensityMatrix nudgeDM <- svToDM <$> genPureSV dim
  let dmw = LA.scale (1 - (weightFraction :+ 0)) (getDensityMatrix dm)
      dmwn = LA.scale (weightFraction :+ 0) nudgeDM
  return $ WeighedDensityMatrix (w, DensityMatrix $ dmw + dmwn)

-- | Calculates values of empirical distribution function at data points.
ecdf :: V.Vector WeighedDensityMatrix -> V.Vector Double
ecdf = V.postscanl' (+) 0 . V.map (fst . getWDM)

-- | /O(log n)/ Given a non-empty sorted vector (typically an empirical cdf
-- evaluated at data points returned by ecdf) and a real number return the
-- (0-based) index of the least element of vector which is greater or equal to
-- the given real number (or the index of the last element, in case there is no
-- element smaller than the argument).
icdf :: V.Vector Double -> Double -> Int
icdf cdf x =
  let tIdx = V.length cdf - 1
      go (lIdx, hIdx) =
        let mIdx =
              truncate $ ((fromIntegral lIdx :: Double) + fromIntegral hIdx) / 2
         in select
              (go (lIdx, mIdx))
              [ (lIdx == hIdx, lIdx)
              , (lIdx + 1 == hIdx, if' (cdf V.! lIdx > x) lIdx hIdx)
              , (cdf V.! mIdx < x, go (mIdx, hIdx))
              ]
   in select (go (0, tIdx)) [(x <= V.head cdf, 0), (x > V.last cdf, tIdx)]

-- | Multinomial resampling of particle vector, which equalises weights of
-- particles.
resampleMultinom :: MWC.GenIO -> Particles -> IO Particles
resampleMultinom gen pts@Particles {..} = do
  us <- MWC.uniformVector gen ptsNumber
  let cdf = ecdf ptsParticles
      idxs = V.map (icdf cdf) us
      w = 1 / fromIntegral ptsNumber
      pointR = over WeighedDensityMatrix (\(_, dm) -> (w, dm))
      resampled = V.map (ptsParticles V.!) idxs
      normed = V.map pointR resampled
  return pts {ptsParticles = normed}

mhmcStep ::
     MWC.GenIO
  -> Dim
  -> Double
  -> [PureStateVector]
  -> Particles
  -> IO (Double, Particles)
mhmcStep gen dim rw ms pts@Particles {..} = do
  let cr wdm wdm' =
        exp
          (fullDataLogLikelihood ms (snd . getWDM $ wdm') -
           fullDataLogLikelihood ms (snd . getWDM $ wdm))
  newParticles <-
    V.mapM (fmap (truncateRank ptsRank) . nudgeParticle dim rw) ptsParticles
  us <- V.replicateM ptsNumber (MWC.uniform gen :: IO Double)
  let crs = V.zipWith cr ptsParticles newParticles
      change = V.zipWith (<=) us crs
      accRate =
        (fromIntegral . V.length . V.filter id) change / fromIntegral ptsNumber
      rwdms = V.zipWith3 if' change newParticles ptsParticles
      final = pts {ptsParticles = rwdms}
  return (accRate, final)

resampleMHMC ::
     ResampleArgs
  -> DensityMatrix
  -> Double
  -> Int
  -> [PureStateVector]
  -> Particles
  -> IO Particles
resampleMHMC ra@ResampleArgs {..} estimate wr iter mts pts = do
  (accRate, resampled) <- mhmcStep argGen argDim wr mts pts
  when (argOut == FullOutput) $
    printf
      "(Weight of new particle: %8.3g, MHMC acceptance rate: %8.3g)\n"
      wr
      accRate
  let (iter', wr') =
        select
          (iter + 1, wr)
          [ (accRate < 1e-2, (0, wr * 0.25))
          , (accRate < 1e-1, (0, wr * 0.5))
          , (iter < argMinIter, (iter + 1, wr))
          , (accRate < 0.33, (0, wr * 0.5))
          ]
  if iter > argMinIter
    then return resampled
    else resampleMHMC ra estimate wr' iter' mts resampled

-- | Arguments for the resampling function.
data ResampleArgs = ResampleArgs
  { argOut :: OutputVerb
  , argGen :: MWC.GenIO
  , argDim :: Dim
  , argMinIter :: MHMCiter
  }

-- | Resample particles. First does one multinomial step that equalises the
-- weights, then performs MHMC iterations adaptively refining the proposal
-- distribution based on acceptance rate. 'argMinIter' iterations are performed
-- for proposal distributions with significant acceptance rates.
resample :: ResampleArgs -> [PureStateVector] -> Particles -> IO Particles
resample ra@ResampleArgs {..} mts pts@Particles {..} = do
  let estimate = getMixedEstimate . V.singleton $ pts
      nudgeW = 0.95
  when (argOut == FullOutput) $ do
    putStrLn ""
    putStrLn $ "resampling rank " ++ show ptsRank
    putStrLn ""
  rm <- resampleMultinom argGen pts
  resampleMHMC ra estimate nudgeW 0 mts rm

max_stars_repo_name: Belinsky-L-V/HABQT

max_issues_repo_name: Belinsky-L-V/HABQT

max_forks_repo_name: Belinsky-L-V/HABQT


text:
{-# language FlexibleInstances      #-}
{-# language KindSignatures         #-}
{-# language FunctionalDependencies #-}
{-# language LambdaCase             #-}
{-# language MultiParamTypeClasses  #-}
{-# language Rank2Types             #-}
{-# language TupleSections          #-}
module Main where

import Control.Monad.Reader
import Data.Complex
import Data.Functor.Const
import Data.Functor.Contravariant
import Data.Functor.Identity
import Data.Monoid                 (First(..))
import Data.Profunctor
import Data.Profunctor.Unsafe      ((#.), (.#))
import Data.Tagged
import Data.Void                   (absurd, vacuous)

type Equality s t a b
  = forall (p :: * -> * -> *) f.
    p a (f b) -> p s (f t)

type Lens      s t a b
  = forall   f. Functor f
  => (a -> f b) -> s -> f t

type Traversal s t a b
  = forall   f. Applicative f
  => (a -> f b) -> s -> f t

type Iso       s t a b
  = forall p f. (Functor f, Profunctor p)
  => p a (f b)  -> p s (f t)

type Prism     s t a b
  = forall p f. (Applicative f, Choice p)
  => p a (f b)  -> p s (f t)

type Getter    s   a
  = forall   f. (Functor f, Contravariant f)
  => (a -> f a) -> s -> f s

type Fold      s   a
  = forall   f. (Applicative f, Contravariant f)
  => (a -> f a) -> s -> f s

type Optic p f s t a b
  = p a (f b)  -> p s (f t)
type Optic' p f t b
  = p b (f b)  -> p t (f t)

type Review t b
  = Tagged b (Identity b) -> Tagged t (Identity t)
-- in other words,
type Review' t b
  = Optic' Tagged Identity t b

-- Traversal:
--
-- We strengthen `Functor` to `Applicative`. Why?
--
-- Clue: note that
--
-- traverse :: Applicative f => (a -> f b) -> t a -> f (t b)
--           ~ Traversable f => Traversal (f a) (f b) a b
--
-- So, instead we can ask why does `traverse` require `Applicative`?
--
-- The instance for lists is informative:
--
-- instance Traversable [] where
--   traverse _ []     = pure []
--   traverse f (x:xs) = (:) <$> f x <*> traverse f xs
--
-- Every lens is a valid Traversal. Witness:

lensToTraversal :: Lens s t a b -> Traversal s t a b
lensToTraversal = id

--
-- Iso:
--
-- We generalize `->` to `(Profunctor p =>) p`. Why?
--
-- One answer is that using `Exchange` (a `Profunctor`) allows us to reverse
-- the `Iso` and extract its parts.
--
-- Prism:
--
-- Here `f` is an `Applicative` and `p` is a `Choice`. Why?
--
-- First of all:
-- * Every `Prism`: Applicative f, Choice p
--   is a valid          v            v
--     `Traversal`: Applicative f
--
--   witness:

prismToTraversal :: Prism s t a b -> Traversal s t a b
prismToTraversal = id

--
--   Here `Prism` adds `Choice p`, meaning TODO
--
-- * Every `Iso`:   Functor f,     Profunctor p
--   is a valid           v             v
--         `Prism`: Applicative f, Choice p
--
--   witness:

isoToPrism :: Iso s t a b -> Prism s t a b
isoToPrism = id

-- Iso's are everything actually:

isoToLens :: Iso s t a b -> Lens s t a b
isoToLens = id

isoToTraversal :: Iso s t a b -> Traversal s t a b
isoToTraversal = id

isoToGetter :: Iso s s a a -> Getter s a
isoToGetter = id

isoToFold :: Iso s s a a -> Fold s a
isoToFold = id

isoToReview :: Iso s s a a -> Review s a
isoToReview = id

-- Except an `Equality`, that is. An `Equality` is an `Iso`:

equalityToIso :: Equality s t a b -> Iso s t a b
equalityToIso = id

-- An Equality is a witness that a ~ s and b ~ t.
--
--
--   Here `Prism` strenthens the constraints:
--   * `Functor` to `Applicative` because it doesn't touch exactly one position
--   * `Profunctor` to `Choice` because TODO
--
-- Getter:
--
-- A getter describes how to retrieve a single value
--
-- > for f to be both Functor and Contravariant implies `f a` doesn't contain
-- > any `a`s at all!
--
-- citation: https://www.reddit.com/r/haskell/comments/5vb6x1/how_do_i_learn_lensinternals/de0uz1v
-- Witness:

coerceGetterF :: (Functor f, Contravariant f) => f a -> f b
coerceGetterF = vacuous . contramap absurd

-- Fold:
--
-- A `Fold` describes how to retrieve multiple values
--
-- Note how conspicuously `Fold` and `Getter` are:
--
-- TODO: understand
-- "A `Getter` is a legal `Fold` that just ignores the supplied `Monoid`
--
-- * Every `Getter`: Functor f,   Contravariant f
--   is a valid          v              v
--           `Fold`: Applicative f, Contravariant f
--
-- Witness:

getterToFold :: Getter s a -> Fold s a
getterToFold = id

-- Like `Getter`, note that the functor used in the `Fold` can't contain any
-- values:

coerceFoldF :: (Applicative f, Contravariant f) => f a -> f b
coerceFoldF = vacuous . contramap absurd

-- A "is a limited form of a `Prism` that can only be used for `re` operations.
-- Witness:

prismToReview :: Prism' t b -> Review t b
prismToReview = id

type Simple f s a  = f s s a a
type Equality' s a = Simple Equality s a
type Lens' s a     = Simple Lens s a
type Iso'  s a     = Simple Iso  s a
type Prism' s a    = Simple Prism s a
type Getting r s a = (a -> Const r a) -> s -> Const r s

iso :: (s -> a) -> (b -> t) -> Iso s t a b
iso sa bt = dimap sa (fmap bt)

-- used to provide access to the two parts of an iso
data Exchange a b s t = Exchange (s -> a) (b -> t)

instance Functor (Exchange a b s) where
  fmap g' (Exchange f g) = Exchange f (g' . g)

instance Profunctor (Exchange a b) where
  dimap f' g' (Exchange f g) = Exchange (f . f') (g' . g)

withIso :: Iso s t a b -> ((s -> a) -> (b -> t) -> r) -> r
withIso ai k = case ai (Exchange id Identity) of
  Exchange sa bt -> k sa (runIdentity #. bt)

-- used to provide access to the two parts of a prism
data Market a b s t = Market (b -> t) (s -> Either t a)

instance Functor (Market a b s) where
  fmap h (Market f g) = Market (h . f) (either (Left . h) Right . g)

instance Profunctor (Market a b) where
  dimap f' g' (Market f g) = Market (g' . f) (left' g' . g . f')

instance Choice (Market a b) where
  left' (Market f g) = Market (Left . f) $ \case
    Left x -> case g x of
      Left y -> Left (Left y)
      Right y -> Right y
    Right c -> Left (Right c)
  -- TODO: implement right'

withPrism :: Prism s t a b -> ((b -> t) -> (s -> Either t a) -> r) -> r
withPrism p k = case p (Market Identity Right) of
  Market bt sa -> k (runIdentity #. bt) (left' runIdentity . sa)

prism :: (b -> t) -> (s -> Either t a) -> Prism s t a b
prism bt seta = dimap seta (either pure (fmap bt)) . right'

matching :: Prism s t a b -> s -> Either t a
matching k = withPrism k $ \_ seta -> seta

_Left :: Prism (Either a c) (Either b c) a b
_Left = prism Left $ either Right (Left . Right)

_Right :: Prism (Either c a) (Either c b) a b
_Right = prism Right $ either (Left . Left) Right

-- TODO: rewrite  / become more comfortable with this section

re :: Review t b -> Getter b t
re p = to (runIdentity #. unTagged #. p .# Tagged .# Identity)

to :: (Profunctor p, Contravariant f) => (s -> a) -> Optic' p f s a
to k = dimap k (contramap k)

reviewSimple :: Prism' s a -> a -> s
reviewSimple r = runIdentity . unTagged . r . Tagged . Identity

review :: MonadReader b m => Optic' Tagged Identity t b -> m t
review r = asks $ runIdentity . unTagged . r . Tagged . Identity

foldMapOf :: Getting r s a -> (a -> r) -> s -> r
foldMapOf l f = getConst #. l (Const #. f)

preview :: Prism' s a -> s -> Maybe a
preview p s = getFirst #. foldMapOf p (First #. Just) $ s

(^?) :: s -> Prism' s a -> Maybe a
(^?) s p = preview p s

--

from :: Iso s t a b -> Iso b a t s
from l = withIso l $ \ sa bt -> iso bt sa

view :: Lens s t a b -> s -> a
view l s = getConst (l Const s)

over :: Lens s t a b -> (a -> b) -> s -> t
over l f = runIdentity . l (Identity . f)

set :: Lens s t a b -> b -> s -> t
set l = over l . const

lens :: (s -> a) -> (s -> b -> t) -> Lens s t a b
lens get' set' f s = set' s <$> f (get' s)

realLens, realLens' :: RealFloat a => Lens' (Complex a) a
realLens f (r :+ i) = fmap (:+ i) (f r)

imagLens, imagLens' :: RealFloat a => Lens' (Complex a) a
imagLens f (r :+ i) = fmap (r :+) (f i)

realLens' = lens (\(r :+ _) -> r) (\(_ :+ i) r -> r :+ i)
imagLens' = lens (\(_ :+ i) -> i) (\(r :+ _) i -> r :+ i)

(<&>) :: Functor f => f a -> (a -> b) -> f b
(<&>) = flip fmap

class Field1 s t a b | s -> a, t -> b, s b -> t, t a -> s where
  _1 :: Lens s t a b

instance Field1 (a, b) (a', b) a a' where
  _1 f ~(a, b) = f a <&> \a' -> (a', b)

instance Field1 (a, b, c) (a', b, c) a a' where
  _1 k ~(a, b, c) = k a <&> \a' -> (a', b, c)

class Field2 s t a b | s -> a, t -> b, s b -> t, t a -> s where
  _2 :: Lens s t a b

instance Field2 (a, b) (a, b') b b' where
  _2 f ~(a, b) = f b <&> \b' -> (a, b')

instance Field2 (a, b, c) (a, b', c) b b' where
  _2 k ~(a, b, c) = k b <&> \b' -> (a, b', c)

main :: IO ()
main = do
  print $ set realLens (1 :: Double) (0 :+ 1)
  print $ view _1 ('x', 'y', 'z')
  print $ view _2 $ set _2 False ('x', 'y', 'z')

max_stars_repo_name: joelburget/mylens

max_issues_repo_name: joelburget/mylens

max_forks_repo_name: joelburget/mylens