text:
{-# LANGUAGE ExtendedDefaultRules #-}
{-# LANGUAGE OverloadedStrings    #-}
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE BangPatterns         #-}

module TestChapter2Bandit (
       testChapter2
      )where

import           Control.Monad
import           Control.Monad.Trans.State (runStateT)
import           Control.Concurrent

import           Data.Configurator
import           Data.Configurator.Types

import           Data.Random
import           Data.Random.Distribution
import           Data.Random.Distribution.Uniform
import           Data.Random.Distribution.Bernoulli
import           Data.Random.Source.IO

import           Data.Text (Text)
import           Data.Time (getCurrentTime)
import           Graphics.Matplotlib

import           Numeric.LinearAlgebra (Vector, Matrix)
import qualified Numeric.LinearAlgebra as LA

import           System.Console.AsciiProgress(Options(..), Stats(..),
                                              displayConsoleRegions, complete,
                                              getProgressStats, def, newProgressBar, tickN)
-- project
import Utils
import Chapter2Bandit

------------------------------------------------------------------------------------------
testChapter2 :: FilePath -> IO ()
testChapter2 configPath = do
  print "Bandit Experiment Starting, will take several minutes "

  (config, _) <- autoReload autoConfig [Required configPath]
  (karm :: Int) <- require config "kArms"
  --
  (bGreedy :: Bool) <- require config "enable.bGreedyExperiment"
  (bUCB :: Bool) <- require config "enable.bUCBExperiment"
  (bGradient :: Bool) <- require config "enable.bGradientExperiment"
  (bCompareTogether :: Bool) <- require config "enable.bCompareTogether"
  (bFigure2_6 :: Bool) <- require config "enable.bFigure2_6"

  -- do the experiments in order
  -- EpsilonGreedy Experiment
  (greedyRewardsCurves, greedyActionCurves) <- if bGreedy then doEpsilonGreedyTests config
                                                          else pure (mp, mp)
  -- UCB Experiment
  (ucbRewardsCurves, ucbActionCurves) <- if bUCB then doUCBTests config else pure (mp, mp)

  -- Gradient Experiment
  (gradientRewardsCurves, gradientActionCurves) <- if bGradient then doGradientTests config
                                                      else pure (mp, mp)

  when bCompareTogether (drawExperimentsCurves config
                              (greedyRewardsCurves % ucbRewardsCurves % gradientRewardsCurves,
                               greedyActionCurves % ucbActionCurves % gradientActionCurves))

  when bFigure2_6 (drawFigure2_6 config)
  pure ()

------------------------------------------------------------------------------------------
-- | doEpsilonGreedyTest Method: compare different epsilons, different initial values,
--   also non-stationary (const step size, weight recent reward more) and stationary environment.
doEpsilonGreedyTests :: Config -> IO (Matplotlib, Matplotlib)
doEpsilonGreedyTests config = do
  startTime <- getCurrentTime
  print (show startTime ++ ": Start Epsilon Greedy Experiment: ")
  (karm :: Int) <- require config "kArms"
  (totalStep :: Int) <- require config "totalStep"
  (totalBandits :: Int) <- require config "totalBandits"
  (ges :: [Double]) <- require config "greedy.epsilons"
  (initValues :: [Double]) <- require config "greedy.initialOptimalValues"
  (stepValues :: [Double]) <- require config "greedy.stepSizes"
  -- run experiments
  let arrParams = zip3 ges initValues stepValues
  !results <- mapM (\ (ge, initValue, stepValue) -> do
                      let !greedyRVar = bernoulli ge
                      !ret <- replicateM totalBandits
                                (oneBanditGo karm totalStep initValue stepValue (EGreedy greedyRVar))
                      pure $ (sum ret) / (fromIntegral totalBandits)
                   )
                  arrParams
  -- draw it
  (ars, oas) <- drawEpsilonGreedy config results
  finishTime <- getCurrentTime
  print $ (show finishTime) ++ ": Finished Epsilon Greedy Experiments"
  pure (ars, oas)

------------------------------------------------------------------------------------------
-- | doUCBTest Method
doUCBTests :: Config -> IO (Matplotlib, Matplotlib)
doUCBTests config = do
  startTime <- getCurrentTime
  print (show startTime ++ ": Start UCB Experiment: ")
  -- read all sample average related params
  (karm :: Int) <- require config "kArms"
  (totalStep :: Int) <- require config "totalStep"
  (totalBandits :: Int) <- require config "totalBandits"
  (sss :: [Double]) <- require config "ucb.stepSizes"
  (eds :: [Double]) <- require config "ucb.exploreDegree"
  -- run experiments
  let arrParams = zip sss eds
  !results <- mapM (\ (stepSize, ed) -> do
                     let policy = UCB ed
                     !ret <- replicateM totalBandits (oneBanditGo karm totalStep 0.0 stepSize policy)
                     pure $ (sum ret) / (fromIntegral totalBandits)
                  ) arrParams
  -- draw it
  (ars, oas) <- drawUCB config results
  finishTime <- getCurrentTime
  print $ (show finishTime) ++ ": Finished UCB Experiments"
  pure (ars, oas)

------------------------------------------------------------------------------------------
-- | doUCBTest Method
doGradientTests :: Config -> IO (Matplotlib, Matplotlib)
doGradientTests config = do
  startTime <- getCurrentTime
  print (show startTime ++ ": Start Gradient Experiment: ")
  -- read all sample average related params
  (karm :: Int) <- require config "kArms"
  (totalStep :: Int) <- require config "totalStep"
  (totalBandits :: Int) <- require config "totalBandits"
  (sss :: [Double]) <- require config "gradient.stepSizes"
  (bBaselines :: [Bool]) <- require config "gradient.baselines"
  -- run experiments
  let arrParams = zip sss bBaselines
  !results <- mapM (\ (stepSize, bBase) -> do
                      let policy = Gradient bBase
                      !ret <- replicateM totalBandits (oneBanditGo karm totalStep 0.0 stepSize policy)
                      pure $ (sum ret) / (fromIntegral totalBandits)
                   ) arrParams
  (ars, oas) <- drawGradient config results
  finishTime <- getCurrentTime
  print $ (show finishTime) ++ ": Finished Gradient Experiments"
  pure (ars, oas)

------------------------------------------------------------------------------------------
oneBanditGo :: Int -> Int -> Double -> Double -> Policy -> IO (Vector Double)
oneBanditGo karm totalStep initValue stepSize policy = do
  trueValues <- generateRandomList karm stdNormal
  let initBandit = mkBandit karm totalStep initValue stepSize trueValues policy
  (averageRewards, bandit) <- runStateT (loopSteps totalStep) initBandit
  -- join two result lists together for convinent
  pure (LA.fromList $ averageRewards ++ (_bestTakes bandit))

------------------------------------------------------------------------------------------
-- drawings
drawFigure2_1 :: Int -> IO ()
drawFigure2_1 karm = do
  kArmRVars::[RVar Double] <- map (flip normal 1.0) <$> generateRandomList karm stdNormal
  kArmDatas <- mapM (generateRandomList 200) kArmRVars
  let figure2_1 = mp % (violinplot kArmDatas @@ [o2 "showmeans" True, o2 "showmedians" True])
                     % xlabel "Action"
                     % ylabel "Reward"
                     % title "Figure 2-1"
                     % grid True
  code figure2_1 >> onscreen figure2_1
  pure ()

drawEpsilonGreedy :: Config -> [Vector Double] -> IO (Matplotlib, Matplotlib)
drawEpsilonGreedy config results = do
  (totalStep :: Int) <- require config "totalStep"
  (ges :: [Double]) <- require config "greedy.epsilons"
  (initValues :: [Double]) <- require config "greedy.initialOptimalValues"

  let rewardCurves = foldl (goPlot 0 totalStep) mp (zip3 ges initValues results)
      bestActions = foldl (goPlot totalStep totalStep) mp (zip3 ges initValues results)
      --subplots @@ [o2 "nrows" 2, o2 "sharey" True]
      -- % setSubplot 0 -- % setSubplot 1
      !figureReward = rewardCurves
                       % title "Epsilon-Greedy Average Reward Comparison"
                       % xlabel "Step"
                       % ylabel "Average Reward"
                       % yticks [0.0, 0.2 .. 1.6]
                       % legend @@ [o2 "fancybox" True,
                                    o2 "shadow" True, o2 "loc" "lower right"]
                       % grid True
      !figureAction = bestActions
                       % title "Epsilon-Greedy Optimal Action Comparison"
                       % xlabel "Step"
                       % ylabel "Optiomal Actions"
                       % yticks [0.0::Double, 0.2 .. 1.0]
                       % legend @@ [o2 "fancybox" True,
                                    o2 "shadow" True, o2 "loc" "lower right"]
                       % grid True
                       % tightLayout
  --
  code figureReward >> code figureAction -- avoid Matplotlib's bug
  onscreen figureReward >> onscreen figureAction
  threadDelay 500000
  pure (rewardCurves, bestActions)
  where
  goPlot :: Int -> Int -> Matplotlib -> (Double, Double, Vector Double) -> Matplotlib
  goPlot startIdx totalLen acc (greedy, initValue, theData) =
    acc % plot [1..totalLen] (LA.subVector startIdx totalLen theData)
               @@ [o2 "label" ("epsilon="++(show greedy)++", Q="++(show initValue))]

drawUCB :: Config -> [Vector Double] -> IO (Matplotlib, Matplotlib)
drawUCB config results = do
  (totalStep :: Int) <- require config "totalStep"
  (eds :: [Double]) <- require config "ucb.exploreDegree"
  let rewardCurves = foldl (goPlot 0 totalStep) mp (zip eds results)
      bestActions = foldl (goPlot totalStep totalStep) mp (zip eds results)
      !figureReward = rewardCurves
                       % title "UCB Average Reward Comparison"
                       % xlabel "Step"
                       % ylabel "Average Reward"
                       % yticks [0.0, 0.2 .. 1.6]
                       % legend @@ [o2 "fancybox" True,
                                    o2 "shadow" True, o2 "loc" "lower right"]
                       % grid True
      !figureAction = bestActions
                       % title "UCB Optimal Action Comparison"
                       % xlabel "Step"
                       % ylabel "Optiomal Actions"
                       % yticks [0.0::Double, 0.2 .. 1.0]
                       % legend @@ [o2 "fancybox" True,
                                    o2 "shadow" True, o2 "loc" "lower right"]
                       % grid True
                       % tightLayout
  --
  code figureReward >> code figureAction -- avoid Matplotlib's bug
  onscreen figureReward >> onscreen figureAction
  threadDelay 500000
  pure (rewardCurves, bestActions)
  where
  goPlot :: Int -> Int -> Matplotlib -> (Double, Vector Double) -> Matplotlib
  goPlot startIdx totalLen acc (exploreDegree, theData) =
    acc % plot [1..totalLen] (LA.subVector startIdx totalLen theData)
               @@ [o2 "label" ("ucb, c=" ++ (show exploreDegree))]

drawGradient :: Config -> [Vector Double] -> IO (Matplotlib, Matplotlib)
drawGradient config results = do
  (totalStep :: Int) <- require config "totalStep"
  (bBaselines :: [Bool]) <- require config "gradient.baselines"
  (sss :: [Double]) <- require config "gradient.stepSizes"
  let rewardCurves = foldl (goPlot 0 totalStep) mp (zip3 sss bBaselines results)
      bestActions = foldl (goPlot totalStep totalStep) mp (zip3 sss bBaselines results)
      !figureReward = rewardCurves
                       % title "Gradient Average Reward Comparison"
                       % xlabel "Step"
                       % ylabel "Average Reward"
                       % yticks [0.0, 0.2 .. 1.6]
                       % legend @@ [o2 "fancybox" True,
                                    o2 "shadow" True, o2 "loc" "lower right"]
                       % grid True
      !figureAction = bestActions
                       % title "Gradient Optimal Action Comparison"
                       % xlabel "Step"
                       % ylabel "Optiomal Actions"
                       % yticks [0.0::Double, 0.2 .. 1.0]
                       % legend @@ [o2 "fancybox" True,
                                    o2 "shadow" True, o2 "loc" "lower right"]
                       % grid True
                       % tightLayout
  --
  code figureReward >> code figureAction -- avoid Matplotlib's bug
  onscreen figureReward >> onscreen figureAction
  threadDelay 500000
  pure (rewardCurves, bestActions)
  where
  goPlot :: Int -> Int -> Matplotlib -> (Double, Bool, Vector Double) -> Matplotlib
  goPlot startIdx totalLen acc (alpha, bBaseline, theData) =
    acc % plot [1..totalLen] (LA.subVector startIdx totalLen theData)
               @@ [o2 "label" ("alpha=" ++ (show alpha) ++ (bBaseline?(",with", ",without")))]

drawExperimentsCurves :: Config -> (Matplotlib, Matplotlib) -> IO ()
drawExperimentsCurves config (rewards, actions) = do
  let !figureReward = rewards
                       % title "Average Reward Comparison"
                       % xlabel "Step"
                       % ylabel "Average Reward"
                       % yticks [0.0, 0.2 .. 1.6]
                       % legend @@ [o2 "fancybox" True,
                                    o2 "shadow" True, o2 "loc" "lower right"]
                       % grid True
      !figureAction = actions
                       % title "Optimal Action Comparison"
                       % xlabel "Step"
                       % ylabel "Optiomal Actions"
                       % yticks [0.0::Double, 0.2 .. 1.0]
                       % legend @@ [o2 "fancybox" True,
                                    o2 "shadow" True, o2 "loc" "lower right"]
                       % grid True
                       % tightLayout
  --
  code figureReward >> code figureAction -- avoid Matplotlib's bug
  onscreen figureReward >> onscreen figureAction
  threadDelay 1000000

------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------
greedyEpsilons = [(-7) .. (-1)]
-- greedy with different initial values and a constant stepSize 0.1
greedyInitValues = [(-2) .. 3]
gradientAlpha  = [(-5) .. 2]
ucbExplore = [(-4) .. 3]

drawFigure2_6 :: Config -> IO ()
drawFigure2_6 config = do
  startTime <- getCurrentTime
  print (show startTime ++ ": Start Figure2_6 Experiment: ")
  -- read all sample average related params
  (karm :: Int) <- require config "kArms"
  (totalStep :: Int) <- require config "totalStep"
  (totalBandits :: Int) <- require config "totalBandits"

  -- 1. different epsilons
  let eGreedyParams = zip3 (map (2.0 ^^) greedyEpsilons) (repeat 0.0) (repeat (-1.0))
  !eGreedyResults <-
    mapM (\ (ge, initValue, stepValue) -> do
            let !greedyRVar = bernoulli ge
            !ret <- replicateM totalBandits
                      (oneBanditGo karm totalStep initValue stepValue (EGreedy greedyRVar))
            pure $ (sum ret) / (fromIntegral totalBandits)
         )  eGreedyParams
  let !eGreedyAverage = calcAverage 0 totalStep eGreedyResults
      !eGreedyAction = calcAverage totalStep totalStep eGreedyResults
      !eGreedyAverageCurve = mp % plot greedyEpsilons eGreedyAverage
                                 @@ [o2 "label" "epsilon-greedy"]
      !eGreedyActionCurve = mp % plot greedyEpsilons eGreedyAction
                                @@ [o2 "label" "epsilon-greedy"]

  -- 2. different e-greedy init values with epsilon = 0.0 and stepSize = 0.1
  let eGreedyInitValParams = zip3 (repeat 0.0) (map (2.0 ^^) greedyInitValues) (repeat 0.1)
  !eGreedyInitValResults <-
    mapM (\ (ge, initValue, stepValue) -> do
            let !greedyRVar = bernoulli ge
            !ret <- replicateM totalBandits
                      (oneBanditGo karm totalStep initValue stepValue (EGreedy greedyRVar))
            pure $ (sum ret) / (fromIntegral totalBandits)
         )  eGreedyInitValParams
  let !eGreedyInitValAverage = calcAverage 0 totalStep eGreedyInitValResults
      !eGreedyInitValActions = calcAverage totalStep totalStep eGreedyInitValResults
      eGreedyInitValAverageCurve = mp % plot greedyInitValues eGreedyInitValAverage
                                             @@ [o2 "label" "greedy-OpticalValue"]
      eGreedyInitValActionCurve = mp % plot greedyInitValues eGreedyInitValActions
                                            @@ [o2 "label" "greedy-OpticalValue"]

  -- 3. ucb with stepSize 0.1
  let ucbParams = zip (repeat 0.1) (map (2.0 ^^) ucbExplore)
  !ucbResults <- mapM (\ (stepSize, ed) -> do
                         let policy = UCB ed
                         !ret <- replicateM totalBandits
                                   (oneBanditGo karm totalStep 0.0 stepSize policy)
                         pure $ (sum ret) / (fromIntegral totalBandits)
                      ) ucbParams
  let !ucbAverage = calcAverage 0 totalStep ucbResults
      !ucbAction = calcAverage totalStep totalStep ucbResults
      ucbAverageCurve = mp % plot ucbExplore ucbAverage @@ [o2 "label" "UCB"]
      ucbActionCurve = mp % plot ucbExplore ucbAction @@ [o2 "label" "UCB"]

  -- 4. gradient method with different alphas and all with baseline
  let gradientParams = zip (map (2.0 ^^) gradientAlpha) (repeat True)
  !gradientResults <- mapM (\ (stepSize, bBase) -> do
                              let policy = Gradient bBase
                              !ret <- replicateM totalBandits
                                        (oneBanditGo karm totalStep 0.0 stepSize policy)
                              pure $ (sum ret) / (fromIntegral totalBandits)
                           ) gradientParams
  let !gradientAverage = calcAverage 0 totalStep gradientResults
      !gradientAction = calcAverage totalStep totalStep gradientResults
      gradientAverageCurve = mp % plot gradientAlpha gradientAverage @@ [o2 "label" "Gradient"]
      gradientActionCurve = mp % plot gradientAlpha gradientAction @@ [o2 "label" "Gradient"]
  -- 5. now draw it
  let !figureReward = eGreedyAverageCurve
                       % gradientAverageCurve
                       % ucbAverageCurve
                       % eGreedyInitValAverageCurve
                       % title "A Parameter study of the various bandit algorithms"
                       % xlabel "Paramter (2^x): e / Alpah / c / Q"
                       % ylabel "Average Reward"
                       % yticks [0.9, 0.1 .. 1.6]
                       % legend @@ [o2 "fancybox" True,
                                    o2 "shadow" True, o2 "loc" "lower right"]
                       % grid True
      !figureAction = eGreedyActionCurve
                       % gradientActionCurve
                       % ucbActionCurve
                       % eGreedyInitValActionCurve
                       % title "A Parameter study of the various bandit algorithms"
                       % xlabel "Paramter (2^x): e / Alpah / c / Q"
                       % ylabel "Optiomal Actions"
                       % yticks [0.2::Double, 0.1 .. 1.0]
                       % legend @@ [o2 "fancybox" True,
                                    o2 "shadow" True, o2 "loc" "lower right"]
                       % grid True
                       % tightLayout
  code figureReward >> code figureAction
  onscreen figureReward >> onscreen figureAction
  threadDelay 500000
  pure ()
  where
  calcAverage startIdx len =
    map ((/ (fromIntegral len)) . sum . LA.toList . LA.subVector startIdx len)

ext: hs

max_stars_repo_name: Xingtao/ReinforceLearningIntro

max_issues_repo_name: Xingtao/ReinforceLearningIntro

max_forks_repo_name: Xingtao/ReinforceLearningIntro


text:
module Q.Risk.VaR where

import Statistics.Autocorrelation
import Data.Csv
import Data.Text (Text)
import qualified Data.ByteString.Lazy as B
import qualified Data.Vector.Storable as VS
import qualified Data.Vector as V
import qualified Numeric.GSL.Statistics as GSL

y :: IO [Double]
y = do
  bytes <- B.readFile "/home/ghais/Downloads/TableF20-1.csv"
  return $ case (decode HasHeader bytes) of
    Left err -> error err
    Right v -> V.toList (V.map fromOnly v)

ar = do
  xs <- demean . VS.fromList <$> y
  let xs2 = VS.map (^2) xs
  return $ GSL.lag1auto  xs2
  where demean xs = let mean = GSL.mean xs
                    in VS.map (\x -> x - mean) xs

max_stars_repo_name: ghais/HQu

max_issues_repo_name: ghais/HQu

max_forks_repo_name: ghais/HQu


text:
module Main where

import           Control.Lens
import           Data.H3.Colour                 (ordinalColours, toCSSColour)
import           Data.H3.Extent                 (Extent, extent, fromTuple,
                                                 resize, toTuple)
import           Data.H3.Scalable               (Scalable (..))
import           Data.H3.Scales                 (IncludeZeroPolicy (..),
                                                 Product (..), continuous)
import           Data.H3.Svg                    (ViewBoxMode (..), renderSvg)
import           Data.H3.Visuals                (ChartVisuals (..),
                                                 FontSize (..),
                                                 LabelOffset (..), Pixel (..),
                                                 Shape (..), Vis (..),
                                                 backgroundShapes, cartesian,
                                                 foregroundShapes, mapColour,
                                                 noGrid)
import           Data.List.NonEmpty             (NonEmpty (..))
import           Data.Semigroup.Foldable        (Foldable1 (..))
import qualified Data.Text                      as Text
import qualified Data.Text.Lazy.IO              as TextIO
import qualified Graphics.Svg                   as SVG
import           Statistics.Distribution        (ContDistr (density),
                                                 Mean (mean), Variance (stdDev))
import qualified Statistics.Distribution.Normal as Distribution.Normal

newtype PlotData = PlotData { getPlots :: [(Double, Double)] }

_PlotData :: Iso' PlotData [(Double, Double)]
_PlotData = iso getPlots PlotData

main :: IO ()
main = do
  let
    plts  = [
      PlotData $ sample 1.0 0.25,
      PlotData $ sample 2.0 0.44]
    dims = (fromTuple (0, 900), fromTuple (500, 0))
    r    = renderSvg (AddViewBox 1.25) (mapColour Text.pack . plotShape plts) dims
  TextIO.writeFile "out.svg" $ SVG.prettyText r

-- | Sample from a normal distribution
sample :: Double -> Double -> [(Double, Double)]
sample theMean theStdDev = vls where
  d = Distribution.Normal.normalDistr theMean theStdDev
  (f, mn, mx) = let m = mean d
                    std = stdDev d in
                    (density d, m - 3 * std, m + 3 * std)
  samples = fmap (\i -> mn + i * delta) [0 .. n] where
    delta = (mx - mn) / n
    n = 100
  vls = fmap (over _2 f) (zip samples samples)

plotShape :: [PlotData] -> (Extent Double, Extent Double) -> Shape String (Pixel Double, Pixel Double)
plotShape ps tgt = result where
  allPoints = foldMap getPlots ps
  result = case allPoints of
    [] -> EmptyShape
    (x:xs) -> over both Pixel <$> AGroup [backgroundShapes s, plotLines, foregroundShapes s] where
      (xExtent, yExtent) = foldMap1 (over both extent) (x :| xs)
      xOpts = continuous xExtent DontIncludeZero
      yOpts = continuous yExtent DontIncludeZero
      tickLength = 5
      ((_, _), (yLw, _)) = over both toTuple tgt
      xyOpts = cartesian
        (fromTuple (yLw, yLw + tickLength))
        (fromTuple (negate tickLength, 0))
        (FontSize "10pt")
        (LabelOffset (-20, 15))
        xOpts
        yOpts
      s = visuals (Vis xyOpts) tgt
      colors = toCSSColour . runIdentity <$> scale (ordinalColours [1 .. length xs]) ()
      plotLines =
        fmap (over both runIdentity . getProduct . combinedScale)
        $ AGroup
        $ mkPlot <$> zip ps [1..]
      mkPlot (PlotData pts, col) = AColouredShape (colors col) $ ALine pts
      combinedScale = scale xyOpts tgt


text:
module Numeric.NLOPTSpec (spec) where

import Numeric.LinearAlgebra
import Numeric.NLOPT
import Test.Hspec

--
-- Next we solve the alkyl problem (see)
--
-- @
--     var x2 := 1.745, >= 0, <= 2;
--     var x3 := 1.2, >= 0, <= 1.6;
--     var x4 := 1.1, >= 0, <= 1.2;
--     var x5 := 3.048, >= 0, <= 5;
--     var x6 := 1.974, >= 0, <= 2;
--     var x7 := 0.893, >= 0.85, <= 0.93;
--     var x8 := 0.928, >= 0.9, <= 0.95;
--     var x9 := 8, >= 3, <= 12;
--     var x10 := 3.6, >= 1.2, <= 4;
--     var x11 := 1.45, >= 1.45, <= 1.62;
--     var x12 := 1, >= 0.99, <= 1.01010101010101;
--     var x13 := 1, >= 0.99, <= 1.01010101010101;
--     var x14 := 1, >= 0.9, <= 1.11111111111111;
--     var x15 := 1, >= 0.99, <= 1.01010101010101;

--     minimize obj:  - 6.3*x5*x8 + 5.04*x2 + 0.35*x3 + x4 + 3.36*x6;

--     subject to

--     e2:  - 0.819672131147541*x2 + x5 - 0.819672131147541*x6 = 0;
--     e3: 0.98*x4 - x7*(0.01*x5*x10 + x4) = 0;
--     e4:  - x2*x9 + 10*x3 + x6 = 0;
--     e5: x5*x12 - x2*(1.12 + 0.13167*x9 - 0.0067*x9*x9) = 0;
--     e6: x8*x13 - 0.01*(1.098*x9 - 0.038*x9*x9) - 0.325*x7 = 0.57425;
--     e7: x10*x14 + 22.2*x11 = 35.82;
--     e8: x11*x15 - 3*x8 = -1.33;
-- @
--
-- Expected solution:
--
-- @
--    x(1) = 3.03582
--    x(2) =  0.95
--    x(3) =  1.7037
--    x(4) =  0.543084
--    x(5) =  0.901319
--    x(6) =  10.4755
--    x(7) =  1.56164
--    x(8) =  1.53535
--    x(9) =  0.99
--    x(10) =  0.99
--    x(11) =  1.11111
--    x(12) =  0.99
--    x(13) =  1.58471
--    x(14) =  2
--    obj =  -1.765
-- @
spec :: Spec
spec = do
{-   describe "Sequential Least-Squares Quadratic Programming (SQSLP)" $ do
    it "solves the Alkyl problem" $ do
      let lb = fromList [0, 0, 0, 0, 0, 0.85, 0.9, 3, 1.2, 1.45, 0.99, 0.99, 0.9, 0.99]
          ub = fromList [2, 1.6, 1.2, 5, 2, 0.93, 0.95, 12, 4, 1.62, 1.01010101010101, 1.01010101010101, 1.11111111111111, 1.01010101010101]
          objfD x =
            ( -6.3 * x ! 3 * x ! 6 + 5.04 * x ! 0 + 0.35 * x ! 1 + x ! 2 + 3.36 * x ! 4,
              fromList [5.04, 0.35, 1, -6.3 * x ! 6, 3.36, 0, -6.3 * x ! 3, 0, 0, 0, 0, 0, 0, 0]
            )
          stop = ObjectiveRelativeTolerance 1e-8 :| []
          eqCons =
            [ EqualityConstraint
                { eqConstraintFunctions =
                    Vector
                      7
                      ( \x _ ->
                          (fromList
                            [ -0.819672131147541 * x ! 0 + x ! 3 - 0.819672131147541 * x ! 4,
                              0.98 * x ! 2 - x ! 5 * (0.01 * x ! 3 * x ! 8 + x ! 2),
                              - x ! 0 * x ! 7 + 10 * x ! 1 + x ! 4,
                              x ! 3 * x ! 10 - x ! 0 * (1.12 + 0.13167 * x ! 7 - 0.0067 * x ! 7 * x ! 7),
                              x ! 6 * x ! 11 - 0.01 * (1.098 * x ! 7 - 0.038 * x ! 7 * x ! 7) - 0.325 * x ! 5 - 0.57425,
                              x ! 8 * x ! 12 + 22.2 * x ! 9 - 35.82,
                              x ! 9 * x ! 13 - 3 * x ! 6 + 1.33
                            ],
                          toDense [
                            ((0,0), -0.819672131147541),
                            ((0,3), 1),
                            ((0,4), -0.819672131147541),
                            ((1,2), 0.98),
                            ((1,5), -(0.01 * x ! 3 * x ! 8 + x ! 2)),
                            ((1,3, -0.01*x ! 5 * x ! 8),
                            ((1,8, -0.01*x ! 5 * x ! 3),
                            ((1,2, -0.01*x ! 5),
                            ((2,0), -x ! 7),
                            ((2,0), -x ! 7),
                            ((2,7), -x ! 0),
                            ((2,1), 10),
                            ((2,4), 1),
                            ((3,3), x ! 10),
                            ((3,10), x ! 3),
                            ((3,0), -(1.12 + 0.13167 * x ! 7 - 0.0067 * x ! 7 * x ! 7)),
                            ((3,7), -x!0*(0.13167 -2*0.0067 * x ! 7)),
                            ((4,))

                                ]
                      ),
                  eqConstraintTolerance = 1e-6
                }
            ]
          subproblem = LocalProblem (fromIntegral $ size lb) stop algorithm
          algorithm = SLSQP objfD [LowerBounds lb, UpperBounds ub] [{- InequalityConstraintsD -}] eqCons
          alkylX0 = [1.75, 1.2, 1.1, 3.048, 1.974, 0.893, 0.982, 8, 3.6, 1.45, 1, 1, 1, 1]
          res = minimizeLocal problem (fromList alkylX0)
      res
        `shouldBe` ( Right
                       ( Solution
                           { solutionCost = -1.765,
                             solutionParams =
                               fromList
                                 [ 3.03582,
                                   0.95,
                                   1.7037,
                                   0.543084,
                                   0.901319,
                                   10.4755,
                                   1.56164,
                                   1.53535,
                                   0.99,
                                   0.99,
                                   1.11111,
                                   0.99,
                                   1.58471,
                                   2
                                 ],
                             solutionResult = FTOL_REACHED
                           }
                       )
                   )
 -}
  describe "Augmented Lagrangian based on SBPLX (Subplex)" $ do
    it "solves the Alkyl problem" $ do
      let lb = fromList [0, 0, 0, 0, 0, 0.85, 0.9, 3, 1.2, 1.45, 0.99, 0.99, 0.9, 0.99]
          ub = fromList [2, 1.6, 1.2, 5, 2, 0.93, 0.95, 12, 4, 1.62, 1.01010101010101, 1.01010101010101, 1.11111111111111, 1.01010101010101]
          objf x = -6.3 * x ! 3 * x ! 6 + 5.04 * x ! 0 + 0.35 * x ! 1 + x ! 2 + 3.36 * x ! 4
          stop = ObjectiveRelativeTolerance 1e-8 :| []
          eqCons =
            [ EqualityConstraint
                { eqConstraintFunctions =
                    Vector
                      7
                      ( \x _ ->
                          fromList
                            [ -0.819672131147541 * x ! 0 + x ! 3 - 0.819672131147541 * x ! 4,
                              0.98 * x ! 2 - x ! 5 * (0.01 * x ! 3 * x ! 8 + x ! 2),
                              - x ! 0 * x ! 7 + 10 * x ! 1 + x ! 4,
                              x ! 3 * x ! 10 - x ! 0 * (1.12 + 0.13167 * x ! 7 - 0.0067 * x ! 7 * x ! 7),
                              x ! 6 * x ! 11 - 0.01 * (1.098 * x ! 7 - 0.038 * x ! 7 * x ! 7) - 0.325 * x ! 5 - 0.57425,
                              x ! 8 * x ! 12 + 22.2 * x ! 9 - 35.82,
                              x ! 9 * x ! 13 - 3 * x ! 6 + 1.33
                            ]
                      ),
                  eqConstraintTolerance = 1e-6
                }
            ]
          problem = AugLagProblem eqCons [] (AUGLAG_EQ_LOCAL subproblem)
          subproblem = LocalProblem (fromIntegral $ size lb) stop algorithm
          algorithm = SBPLX objf [LowerBounds lb, UpperBounds ub] Nothing
          alkylX0 = [1.75, 1.2, 1.1, 3.048, 1.974, 0.893, 0.982, 8, 3.6, 1.45, 1, 1, 1, 1]
          res = minimizeAugLag problem (fromList alkylX0)
      res
        `shouldBe` ( Right
                       ( Solution
                           { solutionCost = -1.765,
                             solutionParams =
                               fromList
                                 [ 3.03582,
                                   0.95,
                                   1.7037,
                                   0.543084,
                                   0.901319,
                                   10.4755,
                                   1.56164,
                                   1.53535,
                                   0.99,
                                   0.99,
                                   1.11111,
                                   0.99,
                                   1.58471,
                                   2
                                 ],
                             solutionResult = FTOL_REACHED
                           }
                       )
                   )


text:
module NN where

import Numeric.LinearAlgebra as LA

rectifiedLinearUnit :: R -> R
rectifiedLinearUnit = max 0

rectifiedLinearUnit' :: R -> R
rectifiedLinearUnit' x = if rectifiedLinearUnit x == 0 then 0 else 1

trainNeuralNetwork :: [((R, R), R)] -> Matrix R -> Matrix R -> Matrix R -> R -> Int -> (Matrix R, Matrix R, Matrix R)
trainNeuralNetwork xys wT b vT learningRate epochs = iterate (step xys learningRate) (wT, b, vT) !! epochs
    where
        step xys learningRate (wT, b, vT) = train xys wT b vT learningRate

train :: [((R, R), R)] -> Matrix R -> Matrix R -> Matrix R -> Double -> (Matrix R, Matrix R, Matrix R)
train xys wT b vT learningRate = train' xys wT b vT learningRate
    where
        train' [] wTs b vTs _ = (wTs, b, vTs)
        train' (((x1,x2), y):xys') wTs b vTs learningRate =
            let (wTs', b', vTs') = trainPoint wTs b vTs learningRate x1 x2 y
            in train' xys' wTs' b' vTs' learningRate

testNeuralNetwork :: [(R, R)] -> Matrix R -> Matrix R -> Matrix R -> R -> [R]
testNeuralNetwork xs wT b vT learningRate = map (calc wT b vT learningRate) xs
    where
        calc wT b vT learningRate (x1, x2) = predict wT b vT learningRate x1 x2

trainPoint :: Matrix R -> Matrix R -> Matrix R -> R -> R -> R -> R -> (Matrix R, Matrix R, Matrix R)
trainPoint wT b vT learningRate x1 x2 y = (wT - scalar learningRate * dw, b - scalar learningRate * db, vT - scalar learningRate * dv)
    where
        x = (2 >< 1) [x1, x2]

        z = wT LA.<> x + b
        a =  cmap rectifiedLinearUnit z

        u =  vT LA.<> a

        du = u - scalar y
        dv = du LA.<> tr a
        da = tr vT LA.<> du
        dz = da * cmap rectifiedLinearUnit' z

        db = dz
        dw = dz LA.<> tr x

predict :: Matrix R -> Matrix R -> Matrix R -> R -> R -> R -> R
predict wT b vT learningRate x1 x2  = u `atIndex` (0,0)
    where
        z = wT LA.<> (2 >< 1) [x1, x2] + b
        a = cmap rectifiedLinearUnit z
        u = vT LA.<> a


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MagicHash             #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE UnboxedTuples         #-}
{-# LANGUAGE ViewPatterns          #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Numeric.Matrix.Internal.Float () where

import qualified Control.Monad.ST        as ST
import qualified Numeric.DataFrame.ST    as ST
import           Numeric.Dimensions      (KnownDim)
import           Numeric.Matrix.Internal
import           Numeric.Matrix.LU
import           Numeric.Scalar.Internal
import           Numeric.Vector.Internal

{-# INLINE mkMat #-}
mkMat ::
  Float -> Float -> Float -> Float ->
  Float -> Float -> Float -> Float ->
  Float -> Float -> Float -> Float ->
  Float -> Float -> Float -> Float ->
  Mat44f
mkMat
  _11 _12 _13 _14
  _21 _22 _23 _24
  _31 _32 _33 _34
  _41 _42 _43 _44
  = ST.runST $ do
    df <- ST.newDataFrame
    ST.writeDataFrameOff df 0  $ scalar _11
    ST.writeDataFrameOff df 1  $ scalar _12
    ST.writeDataFrameOff df 2  $ scalar _13
    ST.writeDataFrameOff df 3  $ scalar _14
    ST.writeDataFrameOff df 4  $ scalar _21
    ST.writeDataFrameOff df 5  $ scalar _22
    ST.writeDataFrameOff df 6  $ scalar _23
    ST.writeDataFrameOff df 7  $ scalar _24
    ST.writeDataFrameOff df 8  $ scalar _31
    ST.writeDataFrameOff df 9  $ scalar _32
    ST.writeDataFrameOff df 10 $ scalar _33
    ST.writeDataFrameOff df 11 $ scalar _34
    ST.writeDataFrameOff df 12 $ scalar _41
    ST.writeDataFrameOff df 13 $ scalar _42
    ST.writeDataFrameOff df 14 $ scalar _43
    ST.writeDataFrameOff df 15 $ scalar _44
    ST.unsafeFreezeDataFrame df

instance HomTransform4 Float where
  {-# INLINE translate4 #-}
  translate4 (unpackV4# -> (# x, y, z, _ #)) = mkMat
    1 0 0 0
    0 1 0 0
    0 0 1 0
    x y z 1

  {-# INLINE translate3 #-}
  translate3 (unpackV3# -> (# x, y, z #)) = mkMat
    1 0 0 0
    0 1 0 0
    0 0 1 0
    x y z 1

  {-# INLINE rotateX #-}
  rotateX a = mkMat
    1 0 0 0
    0 c s 0
    0 n c 0
    0 0 0 1
    where
      c = cos a
      s = sin a
      n = -s

  {-# INLINE rotateY #-}
  rotateY a = mkMat
    c 0 n 0
    0 1 0 0
    s 0 c 0
    0 0 0 1
    where
      c = cos a
      s = sin a
      n = -s

  {-# INLINE rotateZ #-}
  rotateZ a = mkMat
    c s 0 0
    n c 0 0
    0 0 1 0
    0 0 0 1
    where
      c = cos a
      s = sin a
      n = -s

  {-# INLINE rotate #-}
  rotate (unpackV3# -> (# x, y, z #)) a = mkMat
    (c+xxv)  (yxv+zs) (zxv-ys) 0
    (xyv-zs) (c+yyv)  (zyv+xs) 0
    (xzv+ys) (yzv-xs) (c+zzv)  0
     0        0        0       1
    where
      c = cos a
      v = 1 - c -- v for versine
      s = sin a
      xxv = x * x * v
      xyv = x * y * v
      xzv = x * z * v
      yxv = xyv
      yyv = y * y * v
      yzv = y * z * v
      zxv = xzv
      zyv = yzv
      zzv = z * z * v
      xs = x * s
      ys = y * s
      zs = z * s

  {-# INLINE rotateEuler #-}
  rotateEuler x y z = mkMat
    (cy*cz)  (cx*sz+sx*sy*cz) (sx*sz-cx*sy*cz) 0
    (-cy*sz) (cx*cz-sx*sy*sz) (sx*cz+cx*sy*sz) 0
     sy      (-sx*cy)         (cx*cy)          0
     0        0                0               1
    where
      cx = cos x
      sx = sin x
      cy = cos y
      sy = sin y
      cz = cos z
      sz = sin z

  {-# INLINE lookAt #-}
  lookAt up cam foc = mkMat
    xb1 yb1 zb1 0
    xb2 yb2 zb2 0
    xb3 yb3 zb3 0
    tx  ty  tz  1
    where
      (# xb1, xb2, xb3 #) = unpackV3# xb
      (# yb1, yb2, yb3 #) = unpackV3# yb
      (# zb1, zb2, zb3 #) = unpackV3# zb
      zb = normalized $ cam - foc -- Basis vector for "backward", since +Z is behind the camera
      xb = normalized $ up `cross` zb -- Basis vector for "right"
      yb = zb `cross` xb -- Basis vector for "up"
      ncam = -cam
      tx = unScalar $ xb `dot` ncam
      ty = unScalar $ yb `dot` ncam
      tz = unScalar $ zb `dot` ncam

  {-# INLINE perspective #-}
  perspective n f fovy aspect = mkMat
    dpw 0   0   0
    0   dph 0   0
    0   0   a (-1)
    0   0   b   0
    where
      hpd = tan (fovy * 0.5) -- height/distance
      wpd = aspect * hpd; -- width/distance
      dph = recip hpd -- distance/height
      dpw = recip wpd -- distance/width
      nmf = n - f
      a = (n + f) / nmf
      b = 2 * n * f / nmf

  {-# INLINE orthogonal #-}
  orthogonal n f w h = mkMat
    iw 0  0 0
    0  ih 0 0
    0  0  a 0
    0  0  b 1
    where
      ih = 2 / h
      iw = 2 / w
      nmf = n - f
      a = 2 / nmf
      b = (n + f) / nmf

  {-# INLINE toHomPoint #-}
  toHomPoint (unpackV3# -> (# x, y, z #)) = vec4 x y z 1

  {-# INLINE toHomVector #-}
  toHomVector (unpackV3# -> (# x, y, z #)) = vec4 x y z 0

  {-# INLINE fromHom #-}
  fromHom (unpackV4# -> (# x, y, z, w #))
    | w == 0    = vec3 x y z
    | otherwise = vec3 (x/w) (y/w) (z/w)

instance KnownDim n => MatrixInverse Float n where
    inverse = inverseViaLU

instance KnownDim n => MatrixDeterminant Float n where
    det = detViaLU

max_stars_repo_name: achirkin/fasttensor

max_issues_repo_name: achirkin/fasttensor

max_forks_repo_name: achirkin/fasttensor


text:
-- |
-- Module    : Statistics.Matrix.Algorithms
-- Copyright : 2014 Bryan O'Sullivan
-- License   : BSD3
--
-- Useful matrix functions.

module Statistics.Matrix.Algorithms
    (
      qr
    ) where

import Control.Applicative
import Control.Monad.ST (ST, runST)
import Statistics.Matrix (Matrix, column, dimension, for, norm)
import qualified Statistics.Matrix.Mutable as M
import Statistics.Sample.Internal (sum)
import qualified Data.Vector.Unboxed as U
import Prelude hiding (sum, replicate)

-- | /O(r*c)/ Compute the QR decomposition of a matrix.
-- The result returned is the matrices (/q/,/r/).
qr :: Matrix -> (Matrix, Matrix)
qr mat = runST $ do
  let (m,n) = dimension mat
  r <- M.replicate n n 0
  a <- M.thaw mat
  for 0 n $ \j -> do
    cn <- M.immutably a $ \aa -> norm (column aa j)
    M.unsafeWrite r j j cn
    for 0 m $ \i -> M.unsafeModify a i j (/ cn)
    for (j+1) n $ \jj -> do
      p <- innerProduct a j jj
      M.unsafeWrite r j jj p
      for 0 m $ \i -> do
        aij <- M.unsafeRead a i j
        M.unsafeModify a i jj $ subtract (p * aij)
  (,) <$> M.unsafeFreeze a <*> M.unsafeFreeze r

innerProduct :: M.MMatrix s -> Int -> Int -> ST s Double
innerProduct mmat j k = M.immutably mmat $ \mat ->
  sum $ U.zipWith (*) (column mat j) (column mat k)

max_stars_repo_name: runeksvendsen/hs-gauge

max_issues_repo_name: runeksvendsen/hs-gauge

max_forks_repo_name: runeksvendsen/hs-gauge


text:
module Lib
    ( forwardPropogate,
      weightsWithInitializer,
      biasesWithInitializer,
      randomWeights,
      sigmoid,
      runBatch,
      updateParameters,
      getLayer,
      runEpochs,
      ActivationFunction (..),
      Activation (..),
      Parameters,
      Gradients,
      Layer (..)
    ) where

import Numeric.LinearAlgebra

-- |The Layer type represents a layer in an artificial neural network
data Layer = Layer { layerParameters  :: Parameters
                   , activationFunction :: ActivationFunction
                   }

-- |Wraps activation function and derivative
data ActivationFunction = ActivationFunction { f  :: Double -> Double
                                             , f' :: Double -> Double
                                             }
-- |Wraps values from layer activations
data Activation = Activation { -- | Values after matrix multiplication but before activation function is applied
                               preActivation     :: Matrix Double
                               -- | Values after activation function is applied
                             , postActivation    :: Matrix Double
                             }

type Parameters = (Matrix Double, Vector Double)
type Gradients = Parameters

-- |Computes values for each layer
forwardPropogate :: Matrix Double -> [Layer] -> [Activation]
forwardPropogate inputs = scanl activate firstActivation
  where
    firstActivation = Activation {preActivation=(numInputRows><numInputCols) [0..], postActivation=inputs}
    numInputRows = rows inputs
    numInputCols = cols inputs

activate :: Activation -> Layer -> Activation
activate previousActivation layer = Activation {preActivation=layerPreActivation, postActivation=layerPostActivation}
  where
    layerPostActivation = cmap layerActivationFunction layerPreActivation
    layerPreActivation = (previousPostActivation <> layerWeights) `add` fromRows (replicate (rows previousPostActivation) layerBiases)
    previousPostActivation = postActivation previousActivation
    (layerWeights,layerBiases) = layerParameters layer
    layerActivationFunction = (f . activationFunction) layer

-- |Get random weights for a given network structure
randomWeights :: [Int] -> Maybe (IO [Matrix Double])
randomWeights []  = Nothing
randomWeights [x] = Nothing
randomWeights xs  = Just $ traverse (uncurry randn) $ steps xs

-- |Get weights all with initializer value for given network structure
weightsWithInitializer :: Double -> [Int] -> [Matrix Double]
weightsWithInitializer initializer layerSizes = map (\(from,to) -> (from><to) $ repeat initializer) $ steps layerSizes

biasesWithInitializer :: Double -> [Int] -> [Vector Double]
biasesWithInitializer initializer layerSizes = map ((fromList . ($ initializer)) . replicate) $ tail layerSizes

steps :: [a] -> [(a,a)]
steps [] = []
steps [x] = []
steps (x:y:xs) = (x,y) : steps (y:xs)

-- |Sigmoid activation function
sigmoid :: ActivationFunction
sigmoid = ActivationFunction {f=s,f'=s'}
  where
    s x = 1 / (1 + exp (-x))
    s' x = s x * (1 - s x)

-- |Run a batch for training neural network model
runBatch :: Matrix Double -> Matrix Double -> [Layer] -> [Gradients]
runBatch inputs expected initialParameters = runBatch' (outputPreviousActivations:remainingActivations) (reverse initialWeights) (tail reversedActivationDerivatives) outputDelta [(outputWeightGradients, outputBiasGradients)]
  where
    outputBiasGradients = getBiasGradient outputDelta
    outputWeightGradients = getWeightGradients outputDelta outputPreviousActivations
    outputDelta = getOutputDelta expected output outputActivationDerivative
    (output:outputPreviousActivations:remainingActivations) = reverse activations
    outputActivationDerivative = head reversedActivationDerivatives
    activations = forwardPropogate inputs initialParameters
    initialWeights = map (fst . layerParameters) initialParameters
    initialBiases = map (snd . layerParameters) initialParameters
    reversedActivationDerivatives = reverse $ map (f' . activationFunction) initialParameters

runBatch' :: [Activation] -> [Matrix Double] -> [Double -> Double] -> Matrix Double -> [Gradients] -> [Gradients]
runBatch' [input] _ _ _ gradients = gradients
runBatch' reversedActivations reversedWeights activationDerivatives nextDelta gradients = runBatch' (previousActivations:remaining) currentWeights remainingActivationDerivatives currentDelta (gradient:gradients)
  where
    gradient = (getWeightGradients currentDelta previousActivations, getBiasGradient currentDelta)
    currentDelta = getDelta currentActivations nextDelta nextWeights currentActivationDerivative
    (currentActivations:previousActivations:remaining) = reversedActivations
    (nextWeights:currentWeights) = reversedWeights
    (currentActivationDerivative:remainingActivationDerivatives) = activationDerivatives

getDelta :: Activation -> Matrix Double -> Matrix Double -> (Double -> Double) -> Matrix Double
getDelta activation nextDelta nextWeights activationDerivative = (nextDelta <> tr nextWeights) * cmap activationDerivative (preActivation activation)

getBiasGradient :: Matrix Double -> Vector Double
getBiasGradient deltaMatrix = fromList $ map (\column -> sum column / fromIntegral (length column)) columns
  where
    columns = map toList $ toColumns deltaMatrix

getOutputDelta :: Matrix Double -> Activation -> (Double -> Double) -> Matrix Double
getOutputDelta expected output activationDerivative = (postActivation output - expected) * cmap activationDerivative (preActivation output)

getWeightGradients ::  Matrix Double -> Activation -> Matrix Double
getWeightGradients delta previousActivations = tr previousPostActivation <> delta
  where
    previousPostActivation = postActivation previousActivations

updateParameters :: Double -> Double -> [Gradients] -> [Parameters] -> [Parameters]
updateParameters learningRate regularizationTerm = zipWith updateParameters
  where
    updateParameters :: Gradients -> Parameters -> Parameters
    updateParameters (weightGradients,biasGradients) (weights,biases) =
      (
        weights - cmap (*learningRate) (applyRegularization regularizationTerm weightGradients),
        biases - cmap (*learningRate) biasGradients
      )

applyRegularization :: Double -> Matrix Double -> Matrix Double
applyRegularization regularizationTerm weightGradients = weightGradients + cmap (*regularizationTerm) weightGradients

runEpochs :: Int -> Double -> Double -> Matrix Double -> Matrix Double -> [Layer] -> [Parameters]
runEpochs numEpochs learningRate regularizationTerm inputs expected = runEpochs' numEpochs
  where
    runEpochs' 0 layers = map layerParameters layers
    runEpochs' remainingEpochs layers = runEpochs' (remainingEpochs - 1) (runEpoch learningRate regularizationTerm layers inputs expected)

runEpoch :: Double -> Double -> [Layer] -> Matrix Double -> Matrix Double -> [Layer]
runEpoch learningRate regularizationTerm layers inputs expected = updatedLayers
  where
    gradientMatrices = runBatch inputs expected layers
    parameters = map layerParameters layers
    functions = map activationFunction layers
    updatedParameters = updateParameters learningRate regularizationTerm gradientMatrices parameters
    updatedLayers = zipWith getLayer updatedParameters functions

getLayer :: Parameters -> ActivationFunction -> Layer
getLayer parameters activationFunction = Layer{layerParameters=parameters, activationFunction=activationFunction}

max_stars_repo_name: dmurphs/hsnet

max_issues_repo_name: dmurphs/hsnet

max_forks_repo_name: dmurphs/hsnet


text:
-- This module has functions that store pinwheels to disk and read them out to do convolution,
-- for the purpose of handling the case when filter size is greater than RAM size.
module STC.Binary
  ( makePlanBinary
  , computeInitialEigenVectorBinary
  , writeDFTPinwheel
  , convolutionBinary
  ) where

import           Array.UnboxedArray           as AU
import           Control.Monad                as M
import           Control.Monad.IO.Class
import           Control.Monad.Trans.Resource
import           Data.Array.Repa              as R
import           Data.Binary
import           Data.ByteString.Lazy         as BSL
import           Data.ByteString.Lazy         as BS
import           Data.Complex
import           Data.Conduit
import           Control.DeepSeq              (($!!))
import           Data.Conduit.List            as CL
import           Data.List                    as L
import           Data.Vector.Storable         as VS
import           Data.Vector.Unboxed          as VU
import           DFT.Plan
import           Filter.Pinwheel
import           FokkerPlanck.Interpolation
import           FokkerPlanck.Pinwheel        (cutoff)
import           GHC.Float
import           System.IO
import           System.Random
import           Types
import           Utils.Array
import           Utils.Parallel               (ParallelParams (..),
                                               parConduitIO)

makePlanBinary ::
     DFTPlan -> Bool -> FilePath -> Int -> Int -> Int -> Int -> IO DFTPlan
makePlanBinary oldPlan wisdomFlag wisdomFilePath numThetaFreqs numScaleFreqs rows cols = do
  let n = numThetaFreqs * numScaleFreqs * rows * cols
  xs <- M.replicateM n randomIO :: IO [Double]
  ys <- M.replicateM n randomIO :: IO [Double]
  let vecTemp = VS.fromList . L.zipWith (:+) xs $ ys
  when wisdomFlag (importFFTWWisdom wisdomFilePath)
  lock <- getFFTWLock
  plan <-
    fst <$>
    (dft1dGPlan
       lock
       oldPlan
       [cols, rows, numThetaFreqs, numScaleFreqs]
       [0, 1]
       vecTemp >>= \(plan, vec) ->
       idft1dGPlan
         lock
         plan
         [cols, rows, numThetaFreqs, numScaleFreqs]
         [0, 1]
         vec >>= \(plan, vec) ->
         dft1dGPlan
           lock
           plan
           [numThetaFreqs, numScaleFreqs, cols, rows]
           [0, 1]
           vec >>= \(plan, vec) ->
           idft1dGPlan
             lock
             plan
             [numThetaFreqs, numScaleFreqs, cols, rows]
             [0, 1]
             vec)
  exportFFTWWisdom wisdomFilePath
  return plan

computeInitialEigenVectorBinary ::
     Int
  -> Int
  -> [Double]
  -> [Double]
  -> [R2S1RPPoint]
  -> R2Z2Array
computeInitialEigenVectorBinary xLen yLen thetaFreqs scaleFreqs xs =
  let numThetaFreqs = L.length thetaFreqs
      numScaleFreqs = L.length scaleFreqs
      xShift = div xLen 2
      (xMin, xMax) =
        if odd xLen
          then (-xShift, xShift)
          else (-xShift, xShift - 1)
      yShift = div yLen 2
      (yMin, yMax) =
        if odd yLen
          then (-yShift, yShift)
          else (-yShift, yShift - 1)
      vec =
        toUnboxedVector .
        AU.accum (+) 0 ((xMin, yMin), (xMax, yMax)) .
        L.map
          (\(R2S1RPPoint (x, y, _, _)) ->
             ((x, y), 1 / (fromIntegral . L.length $ xs))) $
        xs
  in computeS .
     R.traverse
       (fromUnboxed (Z :. xLen :. yLen) vec)
       (const (Z :. numThetaFreqs :. numScaleFreqs :. xLen :. yLen)) $ \f idx@(Z :. tf :. sf :. i :. j) ->
       if tf == div numThetaFreqs 2 && sf == div numScaleFreqs 2
         then f (Z :. i :. j)
         else 0

{-# INLINE convertDFTPinwheel #-}
convertDFTPinwheel ::
     DFTPlan
  -> R.Array D DIM5 Double
  -> R.Array U DIM1 Double
  -> R.Array U DIM1 Double
  -> Double
  -> Double
  -> (Int, Int)
  -> (Int, Int)
  -> IO BS.ByteString
convertDFTPinwheel plan radialArr thetaFreqs scaleFreqs hollowRadius rMax (rows, cols) (t, s) =
  let (Z :. numThetaFreq :. numScaleFreq :. numTheta0Freq :. numScale0Freq :. _) =
        extent radialArr
      pinwheelArr =
        traverse2
          thetaFreqs
          scaleFreqs
          (\_ _ -> (Z :. numTheta0Freq :. numScale0Freq :. cols :. rows)) $ \ft0 fs0 (Z :. t0 :. s0 :. i :. j) ->
          pinwheelFunc
            (PinwheelHollow0 hollowRadius)
            (thetaFreqs R.! (Z :. t) - ft0 (Z :. t0))
            (scaleFreqs R.! (Z :. s) + fs0 (Z :. s0))
            rMax
            0
            (i - center cols)
            (j - center rows)
      interpolatedPinwheelArr =
        radialCubicInterpolation
          (R.slice radialArr (Z :. t :. s :. All :. All :. All))
          1
          pinwheelArr
  in do xx <-
          fmap
            (encode .
             L.map (\(x :+ y) -> double2Float x :+ double2Float y) . VS.toList) .
          dftExecute
            plan
            (DFTPlanID DFT1DG [cols, rows, numTheta0Freq, numScale0Freq] [0, 1]) .
          VS.convert . toUnboxed . computeS . rotate4D2 . makeFilter2D $
          interpolatedPinwheelArr
        return $!! xx

{-# INLINE writeDFTPinwheelSink #-}
writeDFTPinwheelSink ::
     FilePath -> ConduitT BS.ByteString Void (ResourceT IO) ()
writeDFTPinwheelSink filePath = do
  h <- liftIO $ openBinaryFile filePath WriteMode
  go h
  where
    go h = do
      z <- await
      case z of
        Nothing -> liftIO $ hClose h
        Just bs -> do
          liftIO . BSL.hPut h . encode . BS.length $ bs
          liftIO . BS.hPut h $ bs
          go h

writeDFTPinwheel ::
     ParallelParams
  -> DFTPlan
  -> R.Array D DIM5 Double
  -> Double
  -> Double
  -> [Double]
  -> [Double]
  -> Double
  -> (Int, Int)
  -> FilePath
  -> IO ()
writeDFTPinwheel parallelParams plan radialArr hollowRadius cutoffRadius thetaFreqs scaleFreqs rMax (rows, cols) filePath = do
  let (Z :. numThetaFreq :. numScaleFreq :. numTheta0Freq :. numScale0Freq :. _) =
        extent radialArr
      thetaFreqsArr = fromListUnboxed (Z :. L.length thetaFreqs) thetaFreqs
      scaleFreqsArr = fromListUnboxed (Z :. L.length scaleFreqs) scaleFreqs
  runConduitRes $
    CL.sourceList
      [(t, s) | t <- [0 .. numThetaFreq - 1], s <- [0 .. numScaleFreq - 1]] .|
    parConduitIO
      parallelParams
      (convertDFTPinwheel
         plan
         (cutoff (round cutoffRadius) radialArr)
         thetaFreqsArr
         scaleFreqsArr
         hollowRadius
         rMax
         (rows, cols)) .|
    writeDFTPinwheelSink filePath

{-# INLINE sourceBinary #-}
sourceBinary :: FilePath -> ConduitT () BS.ByteString (ResourceT IO) ()
sourceBinary filePath = do
  h <- liftIO $ openBinaryFile filePath ReadMode
  go h
  where
    go handle = do
      lenBS <- liftIO (BSL.hGet handle 8)
      if BSL.null lenBS
        then liftIO $ hClose handle
        else do
          let len = decode lenBS :: Int
          bs <- liftIO . BS.hGet handle $ len
          yield bs
          go handle

{-# INLINE convolution #-}
convolution ::
     Bool
  -> DFTPlan
  -> DIM4
  -> R.Array U DIM1 Double
  -> VS.Vector (Complex Double)
  -> BS.ByteString
  -> IO (VU.Vector (Complex Double))
convolution True plan dim@(Z :. cols :. rows :. numThetaFreq :. numScaleFreq) thetaFreqs input bs = do
  arr <-
    fmap (fromUnboxed dim . VS.convert) .
    dftExecute
      plan
      (DFTPlanID IDFT1DG [cols, rows, numThetaFreq, numScaleFreq] [0, 1]) .
    VS.zipWith (*) input $
    VS.fromList . L.map (\(x :+ y) -> float2Double x :+ float2Double y) . decode $
    bs
  return $!! toUnboxed . sumS . sumS . R.traverse2 arr thetaFreqs const $
    (\f ft idx@(Z :. _ :. _ :. i :. _) -> f idx * exp (0 :+ ft (Z :. i) * pi))
convolution False plan dim@(Z :. cols :. rows :. numThetaFreq :. numScaleFreq) thetaFreqs input bs = do
  x <-
    fmap (toUnboxed . sumS . sumS . fromUnboxed dim . VS.convert) .
    dftExecute
      plan
      (DFTPlanID IDFT1DG [cols, rows, numThetaFreq, numScaleFreq] [0, 1]) .
    VS.zipWith (*) input .
    VS.fromList . L.map (\(x :+ y) -> float2Double x :+ float2Double y) . decode $
    bs
  return $!! x

convolutionBinary ::
     Bool
  -> ParallelParams
  -> DFTPlan
  -> FilePath
  -> [Double]
  -> R.Array U DIM4 (Complex Double)
  -> IO (R.Array U DIM4 (Complex Double))
convolutionBinary sinkFlag parallelParams plan filterFilePath thetaFreqs inputArr = do
  let thetaFreqsArr = fromListUnboxed (Z :. (L.length thetaFreqs)) thetaFreqs
      (Z :. numThetaFreq :. numScaleFreq :. cols :. rows) = extent inputArr
  inputArrF <-
    dftExecute
      plan
      (DFTPlanID DFT1DG [cols, rows, numThetaFreq, numScaleFreq] [0, 1]) .
    VU.convert . toUnboxed . computeS . rotate4D2 $
    inputArr
  xs <-
    runConduitRes $
    sourceBinary filterFilePath .|
    parConduitIO
      parallelParams
      (convolution
         sinkFlag
         plan
         (Z :. cols :. rows :. numThetaFreq :. numScaleFreq)
         thetaFreqsArr
         inputArrF) .|
    CL.consume
  let arr = fromUnboxed (extent inputArr) . VU.concat $ xs
  return $
    if sinkFlag
      then computeS . R.traverse2 arr thetaFreqsArr const $ \f ft idx@(Z :. i :. _ :. _ :. _) ->
             f idx * exp (0 :+ ft (Z :. i) * pi)
      else arr

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# OPTIONS_GHC -Wno-unused-imports   #-}
{-# OPTIONS_GHC -Wno-unused-top-binds #-}

-- |
-- Module      : AOC.Challenge.Day17
-- License     : BSD3
--
-- Stability   : experimental
-- Portability : non-portable
--
-- Day 17.  See "AOC.Solver" for the types used in this module!
--
-- After completing the challenge, it is recommended to:
--
-- *   Replace "AOC.Prelude" imports to specific modules (with explicit
--     imports) for readability.
-- *   Remove the @-Wno-unused-imports@ and @-Wno-unused-top-binds@
--     pragmas.
-- *   Replace the partial type signatures underscores in the solution
--     types @_ :~> _@ with the actual types of inputs and outputs of the
--     solution.  You can delete the type signatures completely and GHC
--     will recommend what should go in place of the underscores.

module AOC.Challenge.Day17 (
    day17a
  , day17b
  ) where

import           AOC.Prelude

import qualified Data.Graph.Inductive           as G
import qualified Data.IntMap                    as IM
import qualified Data.IntSet                    as IS
import qualified Data.List.NonEmpty             as NE
import qualified Data.List.PointedList          as PL
import qualified Data.List.PointedList.Circular as PLC
import qualified Data.Map                       as M
import qualified Data.OrdPSQ                    as PSQ
import qualified Data.Sequence                  as Seq
import qualified Data.Set                       as S
import qualified Data.Text                      as T
import qualified Data.Vector                    as V
import qualified Linear                         as L
import Safe
import qualified Text.Megaparsec                as P
import Data.Complex
import qualified Text.Megaparsec.Char           as P
import qualified Text.Megaparsec.Char.Lexer     as PP

parseBox :: String -> Maybe (V2 Point)
parseBox = getBounds <=< traverse readMaybe . words . clearOut valid
  where
    valid k = not (isDigit k || k == '-')
    getBounds = \case
      [x1,x2,y1,y2] -> do
        [xMin, xMax] <- pure $ sort [x1,x2]
        [yMin, yMax] <- pure $ sort [y1,y2]
        pure $ V2 (V2 xMin yMin) (V2 xMax yMax)
      _ -> Nothing

-- | Independent tight bounds for each axis
velBounds :: V2 Point -> V2 Point
velBounds (V2 (V2 x1 y1) (V2 x2 _)) = V2 (V2 vx1 y1) (V2 x2 (abs y1))
  where
    vx1 = case quadraticEq (-fromIntegral x1) 0.5 0.5 of
      QReal a b    -> ceiling (max a b)
      QSimul a     -> ceiling a
      QComplex _ _ -> error "invalid box"

data QuadraticSol = QReal    Double Double
                  | QSimul   Double
                  | QComplex (Complex Double) (Complex Double)

quadraticEq :: Double -> Double -> Double -> QuadraticSol
quadraticEq a b c = case compare discr 0 of
    LT -> QComplex (term1 :+ term2) (term1 :+ (-term2))
    EQ -> QSimul term1
    GT -> QReal (term1 + term2) (term1 - term2)
  where
    discr = b*b - 4 * a * c
    term1 = -0.5*b/a
    term2 = 0.5 * sqrt (abs discr) / a

day17a :: V2 Point :~> _
day17a = MkSol
    { sParse = parseBox
    , sShow  = show
    , sSolve = \(bbox@(V2 (V2 _ y1) (V2 x2 _))) ->
        let V2 (V2 vx1 vy1) (V2 vx2 vy2) = velBounds bbox
        in  maximumMay
              [ mx
              | x <- [vx1..vx2]
              , y <- [vy1..vy2]
              , let steps = stepUntilTooDeep x2 y1 (V2 x y)
              , any (inBoundingBox bbox) steps
              , Just mx <- [maximumMay (map (view _y) steps)]
              ]
    }

stepUntilTooDeep :: Int -> Int -> Point -> [Point]
stepUntilTooDeep xmax ymin v0 = map fst . takeWhile p $ iterate simStep (0, v0)
  where
    p (V2 x y, _) = y >= ymin && x <= xmax

simStep :: (Point, Point) -> (Point, Point)
simStep (pos, vel@(V2 vx vy)) = (pos + vel, vel')
  where
    vel' = V2 (vx - signum vx) (vy - 1)

-- target area: x=248..285, y=-85..-56

day17b :: V2 Point :~> _
day17b = MkSol
    { sParse = sParse day17a
    , sShow = show
    , sSolve = \(bbox@(V2 (V2 _ y1) (V2 x2 _))) ->
        let V2 (V2 vx1 vy1) (V2 vx2 vy2) = velBounds bbox
        in  Just $ countTrue (any (inBoundingBox bbox))
              [ stepUntilTooDeep x2 y1 (V2 x y)
              | x <- [vx1..vx2]
              , y <- [vy1..vy2]
              ]
    }


text:
module Fibon.Run.BenchmarkRunner (
    RunResult(..)
  , RunFailure(..)
  , Fibon.Run.BenchmarkRunner.run
)
where

import Control.Concurrent
import Control.Monad
import Control.Exception
import qualified Data.ByteString as B
import Data.Maybe
import Data.Time.Clock
import qualified Data.Vector.Unboxed as Vector
import Fibon.BenchmarkInstance
import Fibon.Result
import Fibon.Run.BenchmarkBundle
import Fibon.Run.Log as Log
import qualified Fibon.Run.SysTools as SysTools
import Statistics.Sample
import System.Directory
import System.FilePath
import System.IO
import System.Process
import Text.Printf

data RunResult =
    Success {runSummary :: RunSummary, runDetails :: [RunDetail]}
  | Failure [RunFailure]
  deriving (Read, Show)

data RunFailure =
    MissingOutput FilePath
  | DiffError     String
  | Timeout
  | ExitError     {exitExpected :: ExitCode, exitActual :: ExitCode}
  deriving (Read, Show)

run :: BenchmarkBundle -> IO RunResult
run bb = do
  let bmk = (bundleName bb)
      pwd = (pathToExeBuildDir bb)
      cmd = (prettyRunCommand bb)
  Log.info $ "Running Benchmark "
  Log.info $ "   BMK: " ++ bmk
  Log.info $ "   PWD: " ++ pwd
  Log.info $ "   CMD: " ++ cmd
  Log.info $ printf "\n@%s|%s|%s" bmk pwd cmd
  runDirect bb

checkResult :: BenchmarkBundle -> ExitCode -> IO (Maybe [RunFailure])
checkResult bb exitCode = do
  outputs <- mapM (checkOutput bb) (output . benchDetails $ bb)
  let results = checkExit bb exitCode : outputs
      errs    = filter isJust results
  case errs of
    [] -> return $ Nothing
    es -> return $ Just (catMaybes es)

checkExit :: BenchmarkBundle -> ExitCode -> Maybe RunFailure
checkExit bb actual = if actual == expected then Nothing else Just ee
  where expected = expectedExit . benchDetails $ bb
        ee       = ExitError {exitExpected = expected, exitActual = actual}

checkOutput :: BenchmarkBundle -> OutputDescription -> IO (Maybe RunFailure)
checkOutput bb (o, Exists) = do
  let f = (destinationToRealFile bb o)
  e <- doesFileExist f
  if e then return   Nothing
       else return $ Just $ MissingOutput ("File "++f++" does not exist")
checkOutput bb (o, Diff diffFile) = do
  e1 <- checkOutput bb (o, Exists)
  e2 <- checkOutput bb (d, Exists)
  e3 <- runDiff f1 f2
  return $ msum [e1, e2, e3]
  where
  d  = OutputFile diffFile
  f1 = (destinationToRealFile bb o)
  f2 = (destinationToRealFile bb d)

runDiff :: FilePath -> FilePath -> IO (Maybe RunFailure)
runDiff f1 f2 = do
  Log.info $ "Diffing files: "++f1++" "++f2
  (r, o, _) <- readProcessWithExitCode (SysTools.diff) [f1, f2] ""
  if r == ExitSuccess then Log.info "No diff error" >>
                           return   Nothing
                      else Log.info "Diff error" >>
                           (return $ Just $ DiffError o)

destinationToRealFile :: BenchmarkBundle -> OutputDestination -> FilePath
destinationToRealFile bb (OutputFile f) = (pathToExeRunDir bb)  </> f
destinationToRealFile bb  Stdout        = (pathToStdoutFile bb)
destinationToRealFile bb  Stderr        = (pathToStderrFile bb)

readExtraStats :: BenchmarkBundle -> IO ExtraStats
readExtraStats bb = do
  let mbStatsFile   = extraStats bb
      statsFile     = fromJust mbStatsFile
      logReadE :: IOException -> IO ExtraStats
      logReadE e =
        Log.warn ("Error reading stats file: "++statsFile++"\n  "++show e)
        >> return B.empty
  case mbStatsFile of
    Nothing -> return B.empty
    Just f  -> do
      handle logReadE $
        bracket (openFile ((pathToExeRunDir bb) </> f) ReadMode)
                (hClose)
                (\h -> B.hGetContents h >>= \s -> B.length s `seq` return s)

type RunStepResult = IO (Either [RunFailure] RunDetail)

runDirect :: BenchmarkBundle -> IO RunResult
runDirect bb = do
  mbDetails <- go count []
  case mbDetails of
    Left e   -> return $ Failure e
    Right ds -> return $ Success (summarize ds) ds
  where
  go 0 ds = return $ Right (reverse ds)
  go n ds = do
    res <- runB bb
    case res of
      Right d -> go (n-1) (d:ds)
      Left e  -> return $ Left e
  runB    = maybe runBenchmarkWithoutTimeout runBenchmarkWithTimeout limit
  limit   = timeout bb
  count   = (iters bb)

summarize :: [RunDetail] -> RunSummary
summarize ds = RunSummary {
      meanTime     = mean times
    , stdDevTime   = stdDev times
    , statsSummary = stats
  }
  where
    times = (Vector.fromList $ map runTime ds)
    stats = case ds of (x:_) -> runStats x; _ -> B.empty

type TimeoutLength = Int
runBenchmarkWithTimeout :: TimeoutLength -> BenchmarkBundle -> RunStepResult
runBenchmarkWithTimeout us bb = do
  resMVar <- newEmptyMVar
  pidMVar <- newEmptyMVar
  tid1 <- forkIO $ (putMVar resMVar . Just) =<< timeBenchmarkExe bb (Just pidMVar)
  _    <- forkIO $ threadDelay us >> putMVar resMVar Nothing
  res <- takeMVar resMVar
  case res of
    Nothing -> do
      Log.info $ "benchmark timed out after "++(show us)++" us"
      -- try to kill the subprocess
      pid <- tryTakeMVar pidMVar
      maybe pass terminateProcess pid
      -- kill the haskell thread
      killThread tid1
      return $ Left [Timeout]
    Just (runDetail, exitCode) -> do
       maybe (Right runDetail) Left `liftM` checkResult bb exitCode

runBenchmarkWithoutTimeout :: BenchmarkBundle -> RunStepResult
runBenchmarkWithoutTimeout bb = do
  (runDetail, exitCode) <- timeBenchmarkExe bb Nothing
  maybe (Right runDetail) Left `liftM` checkResult bb exitCode

timeBenchmarkExe :: BenchmarkBundle            -- benchmark to run
                 -> Maybe (MVar ProcessHandle) -- in case we need to kill it
                 -> IO (RunDetail, ExitCode)
timeBenchmarkExe bb pidMVar = do
  p     <- bundleProcessSpec bb
  start <- getCurrentTime
  (_, _, _, pid) <- createProcess p
  maybe pass (flip putMVar pid) pidMVar
  exit  <- waitForProcess pid
  end   <- getCurrentTime
  mapM_ closeStdIO [std_in  p, std_out p, std_err p]
  stats <- readExtraStats bb
  return $ (RunDetail (realToFrac (diffUTCTime end start)) stats, exit)

closeStdIO :: StdStream -> IO ()
closeStdIO (UseHandle h) = hClose h
closeStdIO _             = return ()

pass :: IO ()
pass = return()

max_stars_repo_name: dmpots/fibon

max_issues_repo_name: dmpots/fibon

max_forks_repo_name: dmpots/fibon


text:
{-# language LambdaCase #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ConstraintKinds #-}
{-# language DerivingVia #-}
{-# language StandaloneDeriving #-}
{-# language ApplicativeDo #-}

-- | This module defines representation of all figures as SVG elements. as well as default style settings.
module Geometry.SVG
  ( Groupable, Group (..)
  , (<+>), group, row, rowSep
  , (<||>), beside, above
  , SVGable (..), SVGContext(..)
  , showSVG
  ) where

import Prelude hiding (writeFile, unwords)
import Graphics.Svg.Core (Attribute (..))
import Graphics.Svg (doctype, svg11_, with, prettyText, (<<-))
import Graphics.Svg (Element, toElement)
import Graphics.Svg.Elements
import Graphics.Svg.Attributes
import Data.Complex
import Data.Monoid
import Data.Maybe
import Data.Text (Text, pack)
import qualified Data.Text.Lazy as LT
import Data.Double.Conversion.Text (toShortest)

import Geometry.Base
import Geometry.Decorations
import Geometry.Point
import Geometry.Circle
import Geometry.Polygon
import Geometry.Plot
import Geometry.Line
import Geometry.Angle
import Geometry.Figures

showt :: Show a => a -> Text
showt = pack . show

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

find :: WithOptions a => (Option -> Maybe b) -> a -> Maybe b
find p d = extractOption p $ defaultOptions d <> options d

extractOption p = getFirst . foldMap (First . p) . getOptions

optFill          = \case {Fill x -> Just x; _ -> Nothing }
optStroke        = \case {Stroke x -> Just x; _ -> Nothing }
optInvisible     = \case {Invisible x -> Just x; _ -> Nothing }
optThickness     = \case {Thickness x -> Just x; _ -> Nothing }
optMultiStroke   = \case {MultiStroke x -> Just x; _ -> Nothing }
optDashing       = \case {Dashing x -> Just x; _ -> Nothing }
optLabelText     = \case {LabelText x -> Just x; _ -> Nothing }
optLabelPosition = \case {LabelPosition x -> Just x; _ -> Nothing }
optLabelOffset   = \case {LabelOffset x -> Just x; _ -> Nothing }
optLabelCorner   = \case {LabelCorner x -> Just x; _ -> Nothing }
optLabelAngle    = \case {LabelAngle x -> Just x; _ -> Nothing }
optSegmentMark   = \case {SegmentMark x -> Just x; _ -> Nothing }

------------------------------------------------------------
{- | Class for objects that could be represented as SVG elements.
Differs from @Graphics.Svg.ToElement@ class in adding
parameter context to a main render function `toSVG`.
-}
class SVGable a where
  toSVG :: a -> SVGContext -> Element
  toSVG _ _ = mempty

  fmtSVG :: a -> Text
  fmtSVG = mempty

instance SVGable Double where
  fmtSVG n = if n ~= 0 then "0" else toShortest n

instance SVGable Cmp where
  fmtSVG = fmtSVG . xy

instance SVGable XY where
  fmtSVG (x, y) = fmtSVG x <> "," <> fmtSVG y <> " "

instance {-# OVERLAPPABLE #-}
  (Foldable t, SVGable a) => SVGable (t a) where
  toSVG = foldMap toSVG

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

instance SVGable a => SVGable (Decorated a) where
  toSVG d = case find optInvisible d of
    Just True -> mempty
    _ -> toSVG (fromDecorated d) . updateOptions (options d)

attributes :: (WithOptions f, Figure f) => f -> SVGContext -> [Attribute]
attributes f ctx = mconcat fmt . getOptions $ opts
  where
    opts = options f <> figureOptions ctx
    fmt = [ attr optStroke Stroke_
          , attr optFill Fill_
          , attr optThickness Stroke_width_
          , attr optDashing Stroke_dasharray_ ]

    attr opt a x = maybeToList $ (\s -> a <<- pack s) <$> findOp opt x
    findOp p = getFirst . foldMap (First . p)

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

instance SVGable Point where
  toSVG p = circle_ . attr <> labelElement p
    where
      attr =
        attributes p <>
        const [ Cx_ <<- fmtSVG (getX p)
              , Cy_ <<- fmtSVG (getY p)
              , R_ <<- "3" ]

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

instance SVGable Label where
  toSVG l = labelElement l . updateOptions (options l)

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

instance SVGable Circle where
  toSVG c = circle_ . attr <> labelElement c
    where
      (x :+ y) = center c
      attr = attributes c <>
             const [ Cx_ <<- fmtSVG x
                   , Cy_ <<- fmtSVG y
                   , R_ <<- fmtSVG (radius c) ]

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

instance SVGable Line where
  toSVG l = do
    bx <- figureBox
    let relabel s = s #: lpos ((s @-> 0.95) - cmp s)
    foldMap (toSVG . relabel) $ clipBy bx l

instance SVGable Ray where
  toSVG r = do
    bx <- figureBox
    let relabel s = s #: lpos ((s @-> 0.95) - cmp s)
    foldMap (toSVG . relabel) $ clipBy bx r

instance SVGable Segment where
  toSVG s = element . attr <> labelElement s
    where
      (a, b) = refPoints s
      element = if isTrivial s then mempty else line_
      attr = attributes s <>
             const [ X1_ <<- fmtSVG (getX a)
                   , Y1_ <<- fmtSVG (getY a)
                   , X2_ <<- fmtSVG (getX b)
                   , Y2_ <<- fmtSVG (getY b) ]

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

instance SVGable Polyline where
  toSVG p | isTrivial p = labelElement p
          | otherwise = polyline_ . attr <> labelElement p
    where
      attr = attributes p <>
             const [ Points_ <<- foldMap fmtSVG (vertices p) ]

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

instance SVGable Polygon where
  toSVG p | isTrivial p = labelElement p
          | otherwise = polygon_ . attr <> labelElement p
    where
      attr = attributes p <>
             const [ Points_ <<- foldMap fmtSVG (vertices p) ]

deriving via Polygon instance SVGable Triangle
deriving via Triangle instance SVGable RightTriangle
deriving via Polygon instance SVGable Rectangle

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

instance (AlmostEq a, Pnt a) => SVGable (Plot a) where
  toSVG = toSVG . asPolyline

instance (AlmostEq a, Pnt a) => SVGable (ClosedPlot a) where
  toSVG = toSVG . asPolyline

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

instance SVGable Angle where
  toSVG an' ctx = toSVG (rays <+> group arc) ctx' <>
                 labelElement an ctx'
    where
      an = reflex an'
      t = extractOption optLabelText (figureOptions ctx)
      labelText = case t of
        Just "#" -> show (angleValue an)
        Just v -> v
        Nothing -> ""
      ctx' = updateOptions (options an <> mkOptions [LabelText labelText]) ctx
      Just ns = extractOption optMultiStroke (figureOptions ctx')
      Angle p s e = an
      rays = mkPolyline [p + cmp s, p, p + cmp e] # scaleAt' p 20
      arc = [ plotManifold an # at' (p + cmp s) # scaleAt' p r
            | i <- [1..ns]
            , let r = 12 + fromIntegral i * 4 ]

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

labelElement :: (WithOptions f, Figure f) => f -> SVGContext -> Element
labelElement ff ctx = case lb of
                   Just s -> text $ toElement s
                   Nothing -> mempty
  where
    opts = figureOptions ctx
    f = Decorated (options ff <> opts, ff)

    lb = find optLabelText f
    loff = fromMaybe 0 $ find optLabelOffset f
    lpos = fromMaybe 0 $ find optLabelPosition f
    (cx, cy) = fromMaybe (0,0) $ find optLabelCorner f

    fontSize = 14
    textWidth = fromIntegral $ length lb
    text = text_ $ [ X_ <<- fmtSVG x
                   , Y_ <<- fmtSVG y
                   , Font_size_ <<- "16"
                   , Font_family_ <<- "CMU Serif"
                   , Font_style_ <<- "italic"
                   , Stroke_ <<- "none"
                   , Fill_ <<- "white"] <> offsetX <> offsetY

    x :+ y = lpos + d
    d = scale (fromIntegral fontSize) loff
    offsetX = case signum cx of
                -1 -> [ Text_anchor_ <<- "start" ]
                0 -> [ Text_anchor_ <<- "middle" ]
                1 -> [ Text_anchor_ <<- "end" ]
    offsetY = case signum cy of
                1 -> [ Dy_ <<- showt (-fontSize `div` 4 -1) ]
                0 -> [ Dy_ <<- showt (fontSize `div` 4 +1) ]
                -1 -> [ Dy_ <<- showt (fontSize - 2) ]

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

-- | Constrain for an object that could be included to a group.
type Groupable a = (SVGable a, Show a, Trans a, Figure a, Eq a)

-- | The group of inhomogeneous `Groupable` objects.
--
-- > let f = G . translate (1,0) . scale 0.7 . rotate 30 <>
-- >         G . translate (1,0) . scale 0.6 . rotate (-45)
-- > in G aSegment # iterate f # take 8 # rotate 90
-- << figs/compose.svg >>
--
data Group where
  EmptyFig :: Group
  G :: Groupable a => a -> Group
  Append :: Group -> Group -> Group

instance Eq Group where
  _ == _ = False

instance Semigroup Group where (<>) = Append
instance Monoid Group where mempty = EmptyFig

infixl 3 <+>
-- | The appending operator for groupable objects.
(<+>) :: (Groupable a, Groupable b) => a -> b -> Group
a <+> b = G a <> G b

instance Trans Group where
  transform _ EmptyFig = EmptyFig
  transform t (G f) = G $ transform t f
  transform t (Append x xs) = Append (transform t x) (transform t xs)

instance Show Group where
  show EmptyFig = mempty
  show (G a) = show a
  show (Append x xs) = show x <> show xs

instance SVGable Group where
  toSVG EmptyFig = mempty
  toSVG (G a) = toSVG a
  toSVG (Append a b) = toSVG a <> toSVG b

instance Figure Group where
  refPoint = left . lower . corner
  isTrivial _ = False
  box EmptyFig = mempty
  box (G f) = box f
  box (Append a b) = box a <> box b

-- | Returns a group of homogeneous list of objects.
--
-- > group $ take 5 $ iterate (rotate 5 . scale 1.1) aTriangle
-- <<figs/groups.svg>>
--
group :: Groupable a => [a] -> Group
group = foldMap G

-- | Places one figure beside the other.
--
-- > aTriangle `beside` aSquare
-- <<figs/beside.svg>>
--
beside :: (Groupable a, Groupable b) => a -> b -> Group
beside f1 f2 = f1 <+> f2 # superpose (refPoint f2) (right . lower . corner $ f1)

-- | The operator form of the `beside` function.
infixl 2 <||>
(<||>) :: (Groupable a, Groupable b) => a -> b -> Group
(<||>) = beside

-- | Places one figure above the other.
--
-- > aTriangle `above` aSquare
-- <<figs/above.svg>>
--
-- > let tr t = t `above` (t `beside` t)
-- > in G aCircle # iterate tr # take 5 # rotate 225 # scaleX 0.6
-- <<figs/serp.svg>>
--
above :: (Groupable a, Groupable b) => a -> b -> Group
above f1 f2 = f2 <+> f1 # superpose (refPoint f2) (left . upper . corner $ f2)

{- | Arranges a homogeneous list of figures in a row

> row [ regularPoly n | n <- [3..7] ]
<<figs/row.svg>>
-}
row :: Groupable f => [f] -> Group
row [] = EmptyFig
row fs = foldl1 beside (G <$> fs)

{- | Arranges a homogeneous list of figures in a row, separated by a given separator.

> rowSep (space 1) [ regularPoly n | n <- [3..7] ]
<<figs/rowSep.svg>>
-}
rowSep :: (Groupable s, Groupable f) => s -> [f] -> Group
rowSep s [] = EmptyFig
rowSep s fs = foldl1 (\g f -> g `beside` s `beside` f) (G <$> fs)
------------------------------------------------------------

-- | The record containing image parameters and options.
data SVGContext = SVGContext
  { imageSize :: Int
  , figureBox :: Rectangle
  , border :: Double
  , figureOptions :: Options }

updateOptions :: Options -> SVGContext -> SVGContext
updateOptions opts ctx = ctx {figureOptions = figureOptions ctx <> opts}

wrapSVG :: Element -> SVGContext -> Element
wrapSVG content ctx =
  doctype <>
  with (svg11_ content)
  [ Version_ <<- "1.1"
  , Width_ <<- showt (border ctx + figureWidth (figureBox ctx))
  , Height_ <<- showt (1.5*border ctx + figureHeight (figureBox ctx))
  , Style_ <<- "background : #444;" ]

-- | Creates a SVG contents for geometric objects.
showSVG :: (Figure a, SVGable a) => Int -> a -> LT.Text
showSVG size obj = prettyText (contents ctx)
  where
    brd = 30
    contents = toSVG obj' >>= wrapSVG
    ctx = SVGContext size fb brd mempty
    obj' = obj
           # superpose p0 (0 :: Cmp)
           # reflect 0
           # scale ((fromIntegral size - brd*3) / ((w `max` h) `min` paperSize))
           # translate' ((brd*1.5, brd*1.5) :: XY)
    fb = boxRectangle obj
         # superpose p0 (0 :: Cmp)
         # reflect 0
         # scale ((fromIntegral size - brd) / ((w `max` h) `min` paperSize))
         # translate' ((brd/2, brd/2) :: XY)
    paperSize = 50
    w = figureWidth obj
    h = figureHeight obj
    p0 = left . upper . corner $ obj

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

max_stars_repo_name: MethaHardworker/geometry

max_issues_repo_name: MethaHardworker/geometry

max_forks_repo_name: MethaHardworker/geometry


text:
#!/usr/bin/env stack
-- stack --install-ghc runghc --resolver lts-11.9 --package backprop-0.2.2.0 --package random --package hmatrix-backprop-0.1.2.1 --package statistics --package lens --package one-liner-instances --package split --package ghc-typelits-natnormalise --package ghc-typelits-knownnat --package hmatrix --package hmatrix-vector-sized --package microlens --package vector-sized --package transformers --package type-combinators -- -O2

-- | Replace the second line with this one to have "./model.hs" open a ghci session
-- stack --install-ghc ghci --resolver lts-11.9 --package backprop-0.2.2.0 --package random --package hmatrix-backprop-0.1.2.1 --package statistics --package lens --package one-liner-instances --package split --package ghc-typelits-natnormalise --package ghc-typelits-knownnat --package hmatrix --package hmatrix-vector-sized --package microlens --package vector-sized --package transformers --package type-combinators

{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE DeriveGeneric                            #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE FlexibleInstances                        #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE LambdaCase                               #-}
{-# LANGUAGE MultiParamTypeClasses                    #-}
{-# LANGUAGE PartialTypeSignatures                    #-}
{-# LANGUAGE PatternSynonyms                          #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeInType                               #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE ViewPatterns                             #-}
{-# OPTIONS_GHC -fno-warn-orphans                     #-}
{-# OPTIONS_GHC -fno-warn-partial-type-signatures     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}
{-# OPTIONS_GHC -fwarn-redundant-constraints          #-}

import           Control.Monad.Trans.State
import           Data.Bifunctor
import           Data.Foldable
import           Data.List
import           Data.List.Split
import           Data.Tuple
import           Data.Type.Option
import           GHC.Generics                          (Generic)
import           GHC.TypeNats
import           Lens.Micro hiding                     ((&))
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop
import           Numeric.LinearAlgebra.Static.Vector
import           Numeric.OneLiner
import           System.Random
import qualified Data.Vector.Storable.Sized            as SVS
import qualified Numeric.LinearAlgebra                 as HU
import qualified Numeric.LinearAlgebra.Static          as H
import qualified Prelude.Backprop                      as B

data a :& b = !a :& !b
  deriving (Show, Generic)
infixr 2 :&

type Model p a b = forall z. Reifies z W
                => BVar z p
                -> BVar z a
                -> BVar z b

linReg :: Model (Double :& Double) Double Double
linReg (a :&& b) x = b * x + a

squaredErrorGrad
    :: (Backprop p, Backprop b, Num b)
    => Model p a b      -- ^ Model
    -> a                -- ^ Observed input
    -> b                -- ^ Observed output
    -> p                -- ^ Parameter guess
    -> p                -- ^ Gradient
squaredErrorGrad f x targ = gradBP $ \p ->
    (f p (auto x) - auto targ) ^ 2

trainModel
    :: (Fractional p, Backprop p, Num b, Backprop b)
    => Model p a b      -- ^ model to train
    -> p                -- ^ initial parameter guess
    -> [(a,b)]          -- ^ list of observations
    -> p                -- ^ updated parameter guess
trainModel f = foldl' $ \p (x,y) -> p - 0.1 * squaredErrorGrad f x y p

trainModelIO
    :: (Fractional p, Backprop p, Num b, Backprop b, Random p)
    => Model p a b      -- ^ model to train
    -> [(a,b)]          -- ^ list of observations
    -> IO p             -- ^ parameter guess
trainModelIO m xs = do
    p0 <- (/ 10) . subtract 0.5 <$> randomIO
    return $ trainModel m p0 xs

testTrainLinReg :: IO (Double :& Double)
testTrainLinReg = trainModelIO linReg (concat (replicate 1000 samps))
  where
    samps = [(1,1),(2,3),(3,5),(4,7),(5,9)]

logistic :: Floating a => a -> a
logistic x = 1 / (1 + exp (-x))

feedForward
    :: (KnownNat i, KnownNat o)
    => Model (L o i :& R o) (R i) (R o)
feedForward (w :&& b) x = w #> x + b

feedForwardLog
    :: (KnownNat i, KnownNat o)
    => Model (L o i :& R o) (R i) (R o)
feedForwardLog (w :&& b) x = logistic (w #> x + b)

testTrainPerceptron :: IO [R 1]
testTrainPerceptron = do
    trained <- trainModelIO feedForwardLog $ take 10000 (cycle samps)
    print trained
    return [ evalBP2 feedForwardLog trained r | (r, _) <- samps ]
  where
    samps = [ (H.vec2 0 0, 0)
            , (H.vec2 1 0, 0)
            , (H.vec2 0 1, 0)
            , (H.vec2 1 1, 1)
            ]

logReg :: Model (Double :& Double) Double Double
logReg ab = logistic . linReg ab

feedForwardLog'
    :: (KnownNat i, KnownNat o)
    => Model (L o i :& R o) (R i) (R o)
feedForwardLog' wb = logistic . feedForward wb

softMax :: (Reifies z W, KnownNat n) => BVar z (R n) -> BVar z (R n)
softMax x = konst (1 / sumElements expx) * expx
  where
    expx = exp x

feedForwardSoftMax
    :: (KnownNat i, KnownNat o)
    => Model (L o i :& R o) (R i) (R o)
feedForwardSoftMax wb = softMax . feedForward wb

(<~)
    :: (Backprop p, Backprop q)
    => Model  p       b c
    -> Model       q  a b
    -> Model (p :& q) a c
(f <~ g) (p :&& q) = f p . g q
infixr 8 <~

testTrainTwoLayer :: IO [R 1]
testTrainTwoLayer = do
    trained <- trainModelIO model (take 10000 (cycle samps))
    print trained
    return [ evalBP2 model trained r | (r, _) <- samps ]
  where
    model :: Model _ (R 2) (R 1)
    model = feedForwardLog' @4 @1 <~ feedForwardLog' @2 @4
    samps = [ (H.vec2 0 0, 0)
            , (H.vec2 1 0, 1)
            , (H.vec2 0 1, 1)
            , (H.vec2 1 1, 0)
            ]

type ModelS p s a b = forall z. Reifies z W
                   => BVar z p
                   -> BVar z a
                   -> BVar z s
                   -> (BVar z b, BVar z s)

ar2 :: ModelS (Double :& (Double :& Double)) Double Double Double
ar2 (c :&& (φ1 :&& φ2)) yLast yLastLast =
    ( c + φ1 * yLast + φ2 * yLastLast, yLast )

fcrnn
    :: (KnownNat i, KnownNat o)
    => ModelS ((L o i :& L o o) :& R o) (R o) (R i) (R o)
fcrnn ((wX :&& wS) :&& b) x s = ( y, logistic y )
  where
    y  = (wX #> x) + (wS #> s) + b

(<*~*)
  :: (Backprop p, Backprop q, Backprop s, Backprop t)
    => ModelS  p        s       b c
    -> ModelS       q        t  a b
    -> ModelS (p :& q) (s :& t) a c
(f <*~* g) (p :&& q) x (s :&& t) = (z, s' :&& t')
  where
    (y, t') = g q x t
    (z, s') = f p y s
infixr 8 <*~*

mapS
    :: (forall z. Reifies z W => BVar z b -> BVar z c)
    -> ModelS p s a b
    -> ModelS p s a c
mapS f g p x = first f . g p x

toS :: Model  p   a b
    -> ModelS p s a b
toS f p x s = (f p x, s)

(<*~)
  :: (Backprop p, Backprop q)
    => Model   p         b c
    -> ModelS       q  s a b
    -> ModelS (p :& q) s a c
(f <*~ g) (p :&& q) x = first (f p) . g q x
infixr 8 <*~

unroll
    :: (Backprop a, Backprop b)
    => ModelS p s  a   b
    -> ModelS p s [a] [b]
unroll f p xs s0 = swap $ B.mapAccumL f' s0 xs
  where
    -- we have to re-arrange the order of arguments and tuple a bit to
    -- match what `mapAccumL` expects
    f' s x = swap (f p x s)

unrollLast
    :: (Backprop a, Backprop b)
    => ModelS p s  a  b
    -> ModelS p s [a] b
unrollLast f = mapS (last . sequenceVar) (unroll f)

unrollLast'
    :: Backprop a
    => ModelS p s  a  b
    -> ModelS p s [a] b
unrollLast' f p xs s0 = foldl' go (undefined, s0) (sequenceVar xs)
  where
    go (_, s) x = f p x s

fixState
    :: s
    -> ModelS p s a b
    -> Model  p   a b
fixState s0 f p x = fst $ f p x (auto s0)

zeroState
    :: Num s
    => ModelS p s a b
    -> Model  p   a b
zeroState = fixState 0

trainState
    :: (Backprop p, Backprop s)
    => ModelS  p    s  a b
    -> Model  (p :& s) a b
trainState f (p :&& s) x = fst $ f p x s

prime
    :: Foldable t
    => ModelS p s a b     -- ^ model
    -> p                  -- ^ parameterization
    -> s                  -- ^ initial state
    -> t a                -- ^ priming input
    -> s                  -- ^ primed state
prime f p = foldl' $ evalBP2 (\s x -> snd $ f (auto p) x s)

feedback
    :: (Backprop a, Backprop s)
    => ModelS p s a a     -- ^ model
    -> p                  -- ^ parameterization
    -> s                  -- ^ initial state
    -> a                  -- ^ initial input
    -> [a]                -- ^ inifinite feedback loop
feedback f p s0 x0 = unfoldr go (x0, s0)
  where
    go (x, s) = Just (x, (y, s'))
      where
        -- 'T2' tuples up a pair of 'BVar's into a 'BVar' of a tuple
        (y, s') = evalBP (uncurry T2 . f (auto p) (auto x)) s

testAR2 :: IO [Double]
testAR2 = do
    trained <- trainModelIO model $ take 10000 samps
    let primed = prime    model0 trained 0      (take 19 series)
        output = feedback model0 trained primed (series !! 19)
    print trained
    return $ take 200 output
  where
    -- sine wave with period 25
    series :: [Double]
    series = [ sin (2 * pi * t / 25) | t <- [0..]              ]
    samps  = [ (init c, last c)      | c <- chunksOf 19 series ]
    model0 :: ModelS _ _ Double Double
    model0 = ar2
    model  :: Model  _   [Double] Double
    model  = zeroState $ unrollLast model0

testRNN :: IO [R 1]
testRNN = do
    trained <- trainModelIO model $ take 100000 samps
    let primed = prime   model0  trained 0      (take 19 series)
        output = feedback model0 trained primed (series !! 19)
    return $ take 200 output
  where
    -- sine wave with period 25
    series :: [H.R 1]
    series = [ H.konst (sin (2 * pi * t / 25)) | t <- [0..]              ]
    samps  = [ (init c, last c)                | c <- chunksOf 19 series ]
    model0 :: ModelS _ _ (R 1) (R 1)
    model0 = feedForward @30 @1
         <*~ mapS logistic (fcrnn @1 @30)
    model  :: Model  _   [R 1] (R 1)
    model  = zeroState $ unrollLast model0

recurrently
    :: (Backprop a, Backprop b)
    => Model  p   (a :& b) b
    -> ModelS p b  a       b
recurrently f p x yLast = (y, y)
  where
    y = f p (x :&& yLast)

recurrentlyWith
    :: (Backprop a, Backprop b)
    => (forall z. Reifies z W => BVar z c -> BVar z b)
    -> Model  p   (a :& b) c
    -> ModelS p b  a       c
recurrentlyWith store f p x yLast = (y, store y)
  where
    y = f p (x :&& yLast)

ffOnSplit
    :: forall i o. (KnownNat i, KnownNat o)
    => Model _ (R i :& R o) (R o)
ffOnSplit p (rI :&& rO) = feedForward p (rI # rO)

fcrnn'
    :: (KnownNat i, KnownNat o)
    => ModelS _ (R o) (R i) (R o)
fcrnn' = recurrentlyWith logistic (\p -> feedForward p . uncurryT (#))

lagged
    :: (KnownNat n, 1 <= n)
    => Model  p       (R (n + 1)) b
    -> ModelS p (R n) Double      b
lagged f p x xLasts = (y, xLasts')
  where
    fullLasts    = xLasts & x
    y            = f p fullLasts
    (_, xLasts') = headTail fullLasts

ar :: (KnownNat n, 1 <= n)
   => ModelS _ (R n) Double Double
ar = lagged (\p -> fst . headTail . feedForward @_ @1 p)

ar2' :: ModelS _ (R 2) Double Double
ar2' = ar @2

main :: IO ()
main = do
    putStrLn "Linear regression"
    print =<< testTrainLinReg

    putStrLn "Single layer perceptron learning AND"
    mapM_ print =<< testTrainPerceptron

    putStrLn "Two-layer ANN learning XOR"
    mapM_ print =<< testTrainTwoLayer

    putStrLn "Sine (AR2)"
    ar2Test <- testAR2
    mapM_ print (take 30 ar2Test)
    writeFile "ar2sin.dat" $ unlines (show <$> ar2Test)

    putStrLn "Sine (RNN)"
    rnnTest <- testRNN
    mapM_ print (take 30 rnnTest)
    writeFile "rnnsin.dat" $ unlines (show . HU.sumElements . H.extract <$> rnnTest)

pattern (:&&) :: (Backprop a, Backprop b, Reifies z W)
              => BVar z a -> BVar z b -> BVar z (a :& b)
pattern x :&& y <- (\xy -> (xy ^^. t1, xy ^^. t2)->(x, y))
  where
    (:&&) = isoVar2 (:&) (\case x :& y -> (x, y))
{-# COMPLETE (:&&) #-}

t1 :: Lens (a :& b) (a' :& b) a a'
t1 f (x :& y) = (:& y) <$> f x

t2 :: Lens (a :& b) (a :& b') b b'
t2 f (x :& y) = (x :&) <$> f y

instance (Num a, Num b) => Num (a :& b) where
    (+)         = gPlus
    (-)         = gMinus
    (*)         = gTimes
    negate      = gNegate
    abs         = gAbs
    signum      = gSignum
    fromInteger = gFromInteger

instance (Fractional a, Fractional b) => Fractional (a :& b) where
    (/) = gDivide
    recip = gRecip
    fromRational = gFromRational

instance (Random a, Random b) => Random (a :& b) where
    random g0 = (x :& y, g2)
      where
        (x, g1) = random g0
        (y, g2) = random g1
    randomR (x0 :& y0, x1 :& y1) g0 = (x :& y, g2)
      where
        (x, g1) = randomR (x0, x1) g0
        (y, g2) = randomR (y0, y1) g1

instance (Backprop a, Backprop b) => Backprop (a :& b)

uncurryT
    :: (Backprop a, Backprop b, Reifies z W)
    => (BVar z a -> BVar z b -> BVar z c)
    -> BVar z (a :& b)
    -> BVar z c
uncurryT f x = f (x ^^. t1) (x ^^. t2)

instance (KnownNat n, KnownNat m) => Random (L n m) where
    random = runState . fmap vecL $ SVS.replicateM (state random)
    randomR (xs,ys) = runState . fmap vecL $ SVS.zipWithM (curry (state . randomR))
        (lVec xs) (lVec ys)

instance (KnownNat n) => Random (R n) where
    random = runState $ vecR <$> SVS.replicateM (state random)
    randomR (xs,ys) = runState . fmap vecR $ SVS.zipWithM (curry (state . randomR))
        (rVec xs) (rVec ys)

max_stars_repo_name: mstksg/inCode

max_issues_repo_name: mstksg/inCode

max_forks_repo_name: mstksg/inCode


text:
-- To run:
--
--   stack build :hardware-examples
--
--   stack build :hardware-trace >& ~/Haskell/concat/hardware/out/o1
--
-- You might also want to use stack's --file-watch flag for automatic recompilation.

{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DeriveFoldable      #-}
{-# LANGUAGE DeriveFunctor       #-}
{-# LANGUAGE DeriveGeneric       #-}
{-# LANGUAGE DeriveTraversable   #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE FlexibleInstances   #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE PatternSynonyms     #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeFamilies        #-}
{-# LANGUAGE TypeOperators       #-}
{-# LANGUAGE TypeSynonymInstances #-}

{-# OPTIONS_GHC -Wall #-}

{-# OPTIONS_GHC -Wno-missing-signatures #-}
{-# OPTIONS_GHC -Wno-unused-imports #-}

{-# OPTIONS_GHC -dsuppress-idinfo #-}
{-# OPTIONS_GHC -dsuppress-uniques #-}
{-# OPTIONS_GHC -dsuppress-module-prefixes #-}

-- {-# OPTIONS_GHC -ddump-simpl #-}

-- {-# OPTIONS_GHC -ddump-rule-rewrites #-}
-- {-# OPTIONS_GHC -fsimpl-tick-factor=250 #-}  -- default 100

-- {-# OPTIONS -fplugin-opt=ConCat.Plugin:trace #-}

{-# OPTIONS_GHC -fno-do-lambda-eta-expansion #-}

module Main where

import Control.Applicative (liftA2,liftA3)
import Data.Complex
import Data.List           (mapAccumL)
import GHC.Float (int2Double)   -- TEMP
import GHC.Generics hiding (S)
-- import ShapedTypes.FFT

import ConCat.Pair
import ConCat.FFT
import ConCat.Misc ((:*),R,sqr,magSqr,Unop,Binop,inNew,inNew2,(:+))
import ConCat.Circuit (GenBuses(..),(:>),Ty(..),Buses(..))
import qualified ConCat.RunCircuit as RC
import ConCat.Syntactic (Syn,render)
import ConCat.AltCat (Ok2,ccc,(:**:)(..),toCcc')
import qualified ConCat.AltCat as A

import ConCat.Rebox () -- necessary for reboxing rules to fire

import ConCat.Hardware.Verilog
import ConCat.Rep
import ConCat.Nat
import ConCat.Shaped

main :: IO ()
main = sequence_
  [ putChar '\n' -- return ()
    -- Unary
  , runVerilog' "neg" $ \ (x :: Int) -> - x  -- Yields bit inversion, not 2's complement!
  , runVerilog' "odd" $ \ (x :: Int) -> x `mod` 2

    -- Binary
  , runVerilog' "adder" $ \ (x :: Int, y :: Int) -> x + y

    -- Conditional
  , runVerilog' "cond" $ \ (p :: Bool, x :: Int, y :: Int) -> if p then x else y

    -- FFT, via functor composition
  -- , runVerilog' "fft_fc_pair" $ \ ( pr :: (UPair (Complex Double)) ) -> fft pr
  , runVerilog' "fft_fc_pair" $ \ ( (x0::(Complex Double)) :# x1 ) -> fft (x0 :# x1)
  -- , runVerilog' "fft_fc_quad" $ \ ( fc :: ( (Pair :.: Pair) (Complex Float) )) -> fft fc
  -- TODO: Understand why only the second alternative, below, works.
  --       (The other two cause a compilation time-out error.)
  -- , runVerilog' "fft_fc_octet" $ \ ( fc :: ( ((Pair :.: Pair) :.: (Pair :.: Pair)) (Complex Float) )) -> fft fc
  , runVerilog' "fft_fc_rb3" $ fft @(RBin N3) @Float
  -- , runVerilog' "fft_fc_rb3" $ fft @((Pair :.: Pair) :.: (Pair :.: Pair)) @Float
  -- , runCirc "fft-rb3" $ ccc $ fft @(RBin N3) @Double
  ]

{--------------------------------------------------------------------
    Testing utilities
--------------------------------------------------------------------}

runVerilog' :: (GenBuses a, GenBuses b) => String -> (a -> b) -> IO ()
runVerilog' _ _ = error "runVerilog' called directly"
{-# NOINLINE runVerilog' #-}
{-# RULES "runVerilog'"
  forall n f. runVerilog' n f = runVerilog n $ toCcc' f #-}

max_stars_repo_name: kenranunderscore/concat

max_issues_repo_name: con-kitty/concat

max_forks_repo_name: con-kitty/concat