text:
{-# LANGUAGE BangPatterns         #-}
{-# LANGUAGE DataKinds            #-}
{-# LANGUAGE DeriveGeneric        #-}
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE GADTs                #-}
{-# LANGUAGE LambdaCase           #-}
{-# LANGUAGE PolyKinds            #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE StandaloneDeriving   #-}
{-# LANGUAGE TemplateHaskell      #-}
{-# LANGUAGE TypeApplications     #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE ViewPatterns         #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

import           Control.DeepSeq
import           Criterion.Main
import           Criterion.Types
import           Data.Char
import           Data.Functor.Identity
import           Data.Time
import           GHC.Generics                 (Generic)
import           GHC.TypeLits
import           Lens.Micro
import           Lens.Micro.TH
import           Numeric.Backprop
import           Numeric.Backprop.Class
import           Numeric.LinearAlgebra.Static
import           System.Directory
import qualified Data.Vector                  as V
import qualified Numeric.LinearAlgebra        as HM
import qualified System.Random.MWC            as MWC

type family HKD f a where
    HKD Identity a = a
    HKD f        a = f a

data Layer' i o f =
    Layer { _lWeights :: !(HKD f (L o i))
          , _lBiases  :: !(HKD f (R o))
          }
  deriving (Generic)

type Layer i o = Layer' i o Identity

deriving instance (KnownNat i, KnownNat o) => Show (Layer i o)
instance NFData (Layer i o)

makeLenses ''Layer'

data Network' i h1 h2 o f =
    Net { _nLayer1 :: !(HKD f (Layer i  h1))
        , _nLayer2 :: !(HKD f (Layer h1 h2))
        , _nLayer3 :: !(HKD f (Layer h2 o ))
        }
  deriving (Generic)

type Network i h1 h2 o = Network' i h1 h2 o Identity

deriving instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Show (Network i h1 h2 o)
instance NFData (Network i h1 h2 o)

makeLenses ''Network'

main :: IO ()
main = do
    g     <- MWC.initialize
           . V.fromList
           . map (fromIntegral . ord)
           $ "hello world"
    test0 <- MWC.uniformR @(R 784, R 10) ((0,0),(1,1)) g
    net0  <- MWC.uniformR @(Network 784 300 100 10) (-0.5, 0.5) g
    t     <- getZonedTime
    let tstr = formatTime defaultTimeLocale "%Y%m%d-%H%M%S" t
    createDirectoryIfMissing True "bench-results"
    defaultMainWith defaultConfig
          { reportFile = Just $ "bench-results/mnist-bench_" ++ tstr ++ ".html"
          , timeLimit  = 10
          } [
        bgroup "gradient"
          [ let runTest x y     = gradNetManual x y net0
            in  bench "manual"  $ nf (uncurry runTest) test0
          , let runTest x y     = gradBP (netErr x y) net0
            in  bench "bp-lens" $ nf (uncurry runTest) test0
          , let runTest x y     = gradBP (netErrHKD x y) net0
            in  bench "bp-hkd"  $ nf (uncurry runTest) test0
          , let runTest x y     = gradBP (\n' -> netErrHybrid n' y x) net0
            in  bench "hybrid"  $ nf (uncurry runTest) test0
          ]
      , bgroup "descent"
          [ let runTest x y     = trainStepManual 0.02 x y net0
            in  bench "manual"  $ nf (uncurry runTest) test0
          , let runTest x y     = trainStep 0.02 x y net0
            in  bench "bp-lens" $ nf (uncurry runTest) test0
          , let runTest x y     = trainStepHKD 0.02 x y net0
            in  bench "bp-hkd"  $ nf (uncurry runTest) test0
          , let runTest x y     = trainStepHybrid 0.02 x y net0
            in  bench "hybrid"  $ nf (uncurry runTest) test0
          ]
      , bgroup "run"
          [ let runTest         = runNetManual net0
            in  bench "manual"  $ nf runTest (fst test0)
          , let runTest x       = evalBP (`runNetwork` x) net0
            in  bench "bp-lens" $ nf runTest (fst test0)
          , let runTest x       = evalBP (`runNetworkHKD` x) net0
            in  bench "bp-hkd"  $ nf runTest (fst test0)
          , let runTest x       = evalBP (`runNetHybrid` x) net0
            in  bench "hybrid"  $ nf runTest (fst test0)
          ]
      ]

-- ------------------------------
-- - "Backprop" Lens Mode       -
-- ------------------------------

runLayer
    :: (KnownNat i, KnownNat o, Reifies s W)
    => BVar s (Layer i o)
    -> BVar s (R i)
    -> BVar s (R o)
runLayer l x = (l ^^. lWeights) #>! x + (l ^^. lBiases)
{-# INLINE runLayer #-}

softMax :: (KnownNat n, Reifies s W) => BVar s (R n) -> BVar s (R n)
softMax x = konst' (1 / sumElements' expx) * expx
  where
    expx = exp x
{-# INLINE softMax #-}

runNetwork
    :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o, Reifies s W)
    => BVar s (Network i h1 h2 o)
    -> R i
    -> BVar s (R o)
runNetwork n = softMax
             . runLayer (n ^^. nLayer3)
             . logistic
             . runLayer (n ^^. nLayer2)
             . logistic
             . runLayer (n ^^. nLayer1)
             . auto
{-# INLINE runNetwork #-}

crossEntropy
    :: (KnownNat n, Reifies s W)
    => R n
    -> BVar s (R n)
    -> BVar s Double
crossEntropy t r = negate $ log r <.>! auto t
{-# INLINE crossEntropy #-}

netErr
    :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o, Reifies s W)
    => R i
    -> R o
    -> BVar s (Network i h1 h2 o)
    -> BVar s Double
netErr x t n = crossEntropy t (runNetwork n x)
{-# INLINE netErr #-}

trainStep
    :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o)
    => Double
    -> R i
    -> R o
    -> Network i h1 h2 o
    -> Network i h1 h2 o
trainStep r !x !t !n = n - realToFrac r * gradBP (netErr x t) n
{-# INLINE trainStep #-}

-- ------------------------------
-- - "Backprop" HKD Mode        -
-- ------------------------------

runLayerHKD
    :: (KnownNat i, KnownNat o, Reifies s W)
    => BVar s (Layer i o)
    -> BVar s (R i)
    -> BVar s (R o)
runLayerHKD (splitBV->Layer w b) x = w #>! x + b
{-# INLINE runLayerHKD #-}

runNetworkHKD
    :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o, Reifies s W)
    => BVar s (Network i h1 h2 o)
    -> R i
    -> BVar s (R o)
runNetworkHKD (splitBV->Net l1 l2 l3) = softMax
                                      . runLayerHKD l3
                                      . logistic
                                      . runLayerHKD l2
                                      . logistic
                                      . runLayerHKD l1
                                      . auto
{-# INLINE runNetworkHKD #-}

netErrHKD
    :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o, Reifies s W)
    => R i
    -> R o
    -> BVar s (Network i h1 h2 o)
    -> BVar s Double
netErrHKD x t n = crossEntropy t (runNetworkHKD n x)
{-# INLINE netErrHKD #-}

trainStepHKD
    :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o)
    => Double
    -> R i
    -> R o
    -> Network i h1 h2 o
    -> Network i h1 h2 o
trainStepHKD r !x !t !n = n - realToFrac r * gradBP (netErrHKD x t) n
{-# INLINE trainStepHKD #-}

-- ------------------------------
-- - "Manual" Mode              -
-- ------------------------------

runLayerManual
    :: (KnownNat i, KnownNat o)
    => Layer i o
    -> R i
    -> R o
runLayerManual l x = (l ^. lWeights) #> x + (l ^. lBiases)
{-# INLINE runLayerManual #-}

softMaxManual :: KnownNat n => R n -> R n
softMaxManual x = konst (1 / sumElements expx) * expx
  where
    expx = exp x
{-# INLINE softMaxManual #-}

runNetManual
    :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o)
    => Network i h1 h2 o
    -> R i
    -> R o
runNetManual n = softMaxManual
               . runLayerManual (n ^. nLayer3)
               . logistic
               . runLayerManual (n ^. nLayer2)
               . logistic
               . runLayerManual (n ^. nLayer1)
{-# INLINE runNetManual #-}

gradNetManual
    :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o)
    => R i
    -> R o
    -> Network i h1 h2 o
    -> Network i h1 h2 o
gradNetManual x t (Net (Layer w1 b1) (Layer w2 b2) (Layer w3 b3)) =
    let y1 = w1 #> x
        z1 = y1 + b1
        x2 = logistic z1
        y2 = w2 #> x2
        z2 = y2 + b2
        x3 = logistic z2
        y3 = w3 #> x3
        z3 = y3 + b3
        o0 = exp z3
        o1 = HM.sumElements (extract o0)
        o2 = o0 / konst o1
        -- o3 = - (log o2 <.> t)
        dEdO3 = 1
        dEdO2 = dEdO3 * (- t / o2)
        dEdO1 = - (dEdO2 <.> o0) / (o1 ** 2)
        dEdO0 = konst dEdO1 + dEdO2 / konst o1
        dEdZ3 = dEdO0 * o0
        dEdY3 = dEdZ3
        dEdX3 = tr w3 #> dEdY3
        dEdZ2 = dEdX3 * (x3 * (1 - x3))
        dEdY2 = dEdZ2
        dEdX2 = tr w2 #> dEdY2
        dEdZ1 = dEdX2 * (x2 * (1 - x2))
        dEdY1 = dEdZ1
        dEdB3 = dEdZ3
        dEdW3 = dEdY3 `outer` x3
        dEdB2 = dEdZ2
        dEdW2 = dEdY2 `outer` x2
        dEdB1 = dEdZ1
        dEdW1 = dEdY1 `outer` x
    in  Net (Layer dEdW1 dEdB1) (Layer dEdW2 dEdB2) (Layer dEdW3 dEdB3)
{-# INLINE gradNetManual #-}

trainStepManual
    :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o)
    => Double
    -> R i
    -> R o
    -> Network i h1 h2 o
    -> Network i h1 h2 o
trainStepManual r !x !t !n =
    let gN = gradNetManual x t n
    in  n - (realToFrac r * gN)

-- ------------------------------
-- - "Hybrid" Mode              -
-- ------------------------------

layerOp :: (KnownNat i, KnownNat o) => Op '[Layer i o, R i] (R o)
layerOp = op2 $ \(Layer w b) x ->
    ( w #> x + b
    , \g -> (Layer (g `outer` x) g, tr w #> g)
    )
{-# INLINE layerOp #-}

logisticOp
    :: Floating a
    => Op '[a] a
logisticOp = op1 $ \x ->
    let lx = logistic x
    in  (lx, \g -> lx * (1 - lx) * g)
{-# INLINE logisticOp #-}

softMaxOp
    :: KnownNat n
    => Op '[R n] (R n)
softMaxOp = op1 $ \x ->
    let expx   = exp x
        tot    = sumElements expx
        invtot = 1 / tot
        res    = konst invtot * expx
    in  ( res
        , \g -> res - konst (invtot ** 2) * exp (2 * x) * g
        )
{-# INLINE softMaxOp #-}

softMaxCrossEntropyOp
    :: KnownNat n
    => R n
    -> Op '[R n] Double
softMaxCrossEntropyOp targ = op1 $ \x ->
    let expx   = exp x
        sm     = konst (1 / sumElements expx) * expx
        ce     = negate $ log sm <.> targ
    in  ( ce
        , \g -> (sm - targ) * konst g
        )
{-# INLINE softMaxCrossEntropyOp #-}

runNetHybrid
    :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o, Reifies s W)
    => BVar s (Network i h1 h2 o)
    -> R i
    -> BVar s (R o)
runNetHybrid n = liftOp1 softMaxOp
               . liftOp2 layerOp (n ^^. nLayer3)
               . liftOp1 logisticOp
               . liftOp2 layerOp (n ^^. nLayer2)
               . liftOp1 logisticOp
               . liftOp2 layerOp (n ^^. nLayer1)
               . auto
{-# INLINE runNetHybrid #-}

netErrHybrid
    :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o, Reifies s W)
    => BVar s (Network i h1 h2 o)
    -> R o
    -> R i
    -> BVar s Double
netErrHybrid n t = liftOp1 (softMaxCrossEntropyOp t)
                 . liftOp2 layerOp (n ^^. nLayer3)
                 . liftOp1 logisticOp
                 . liftOp2 layerOp (n ^^. nLayer2)
                 . liftOp1 logisticOp
                 . liftOp2 layerOp (n ^^. nLayer1)
                 . auto
{-# INLINE netErrHybrid #-}

trainStepHybrid
    :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o)
    => Double
    -> R i
    -> R o
    -> Network i h1 h2 o
    -> Network i h1 h2 o
trainStepHybrid r !x !t !n =
    let gN = gradBP (\n' -> netErrHybrid n' t x) n
    in  n - (realToFrac r * gN)
{-# INLINE trainStepHybrid #-}

-- ------------------------------
-- - Operations                 -
-- ------------------------------

infixr 8 #>!
(#>!)
    :: (KnownNat m, KnownNat n, Reifies s W)
    => BVar s (L m n)
    -> BVar s (R n)
    -> BVar s (R m)
(#>!) = liftOp2 . op2 $ \m v ->
  ( m #> v, \g -> (g `outer` v, tr m #> g) )
{-# INLINE (#>!) #-}

infixr 8 <.>!
(<.>!)
    :: (KnownNat n, Reifies s W)
    => BVar s (R n)
    -> BVar s (R n)
    -> BVar s Double
(<.>!) = liftOp2 . op2 $ \x y ->
  ( x <.> y, \g -> (konst g * y, x * konst g)
  )
{-# INLINE (<.>!) #-}

konst'
    :: (KnownNat n, Reifies s W)
    => BVar s Double
    -> BVar s (R n)
konst' = liftOp1 . op1 $ \c -> (konst c, HM.sumElements . extract)
{-# INLINE konst' #-}

sumElements :: KnownNat n => R n -> Double
sumElements = HM.sumElements . extract
{-# INLINE sumElements #-}

sumElements'
    :: (KnownNat n, Reifies s W)
    => BVar s (R n)
    -> BVar s Double
sumElements' = liftOp1 . op1 $ \x -> (sumElements x, konst)
{-# INLINE sumElements' #-}

logistic :: Floating a => a -> a
logistic x = 1 / (1 + exp (-x))
{-# INLINE logistic #-}

-- ------------------------------
-- - Instances                  -
-- ------------------------------

instance (KnownNat i, KnownNat o) => Num (Layer i o) where
    Layer w1 b1 + Layer w2 b2 = Layer (w1 + w2) (b1 + b2)
    Layer w1 b1 - Layer w2 b2 = Layer (w1 - w2) (b1 - b2)
    Layer w1 b1 * Layer w2 b2 = Layer (w1 * w2) (b1 * b2)
    abs    (Layer w b)        = Layer (abs    w) (abs    b)
    signum (Layer w b)        = Layer (signum w) (signum b)
    negate (Layer w b)        = Layer (negate w) (negate b)
    fromInteger x             = Layer (fromInteger x) (fromInteger x)

instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Num (Network i h1 h2 o) where
    Net a b c + Net d e f = Net (a + d) (b + e) (c + f)
    Net a b c - Net d e f = Net (a - d) (b - e) (c - f)
    Net a b c * Net d e f = Net (a * d) (b * e) (c * f)
    abs    (Net a b c)    = Net (abs    a) (abs    b) (abs    c)
    signum (Net a b c)    = Net (signum a) (signum b) (signum c)
    negate (Net a b c)    = Net (negate a) (negate b) (negate c)
    fromInteger x         = Net (fromInteger x) (fromInteger x) (fromInteger x)

instance (KnownNat i, KnownNat o) => Fractional (Layer i o) where
    Layer w1 b1 / Layer w2 b2 = Layer (w1 / w2) (b1 / b2)
    recip (Layer w b)         = Layer (recip w) (recip b)
    fromRational x            = Layer (fromRational x) (fromRational x)

instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Fractional (Network i h1 h2 o) where
    Net a b c / Net d e f = Net (a / d) (b / e) (c / f)
    recip (Net a b c)     = Net (recip a) (recip b) (recip c)
    fromRational x        = Net (fromRational x) (fromRational x) (fromRational x)

instance KnownNat n => MWC.Variate (R n) where
    uniform g = randomVector <$> MWC.uniform g <*> pure Uniform
    uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g

instance (KnownNat m, KnownNat n) => MWC.Variate (L m n) where
    uniform g = uniformSample <$> MWC.uniform g <*> pure 0 <*> pure 1
    uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g

instance (KnownNat i, KnownNat o) => MWC.Variate (Layer i o) where
    uniform g = Layer <$> MWC.uniform g <*> MWC.uniform g
    uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g

instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => MWC.Variate (Network i h1 h2 o) where
    uniform g = Net <$> MWC.uniform g <*> MWC.uniform g <*> MWC.uniform g
    uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g

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

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

instance (KnownNat i, KnownNat o) => Backprop (Layer i o)
instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Backprop (Network i h1 h2 o)


text:
module Cas.Instances.Other where

import PreludeCustom
import Cas.Internal.Interface
import Cas.Internal.Normalization
import Data.Complex.Generic

instance Enum X where
    fromEnum = xI
    toEnum = X

instance Fractional T where
    recip = normalize . mInv
    fromRational = TF . fromRational

instance Num T where
    negate      = normalize . aInv
    (+)         = (normalize .) .  addBin
    (*)         = (normalize .) .  mulBin
    fromInteger = TF . fromInteger

instance Floating T where
    (**) = (normalize .) . pow

-- instance Fractional T where
--     recip = mInv
--     fromRational = TF . fromRational

-- instance Num T where
--     negate      =  aInv
--     (+)         =  addBin
--     (*)         =  mulBin
--     fromInteger = TF . fromInteger

-- instance Floating T where
--     (**) = pow

ext: hs

max_stars_repo_name: vcanadi/cas

max_issues_repo_name: vcanadi/cas

max_forks_repo_name: vcanadi/cas


text:
{-# LANGUAGE OverloadedStrings #-}
module VisualizeFunction where
import Control.Monad
import Data.Binary.Put
import Numeric.LinearAlgebra.HMatrix
import System.IO
import qualified Data.ByteString.Lazy as L
import qualified Data.Vector.Storable as V

-- derived from description at https://en.wikipedia.org/wiki/BMP_file_format
putBMPHeader headerLength imageLength = do
    putLazyByteString "BM" -- magic identifier
    putWord32le (headerLength + imageLength) -- total length of file (bytes)
    putLazyByteString "HELO" -- apparently a reserved field which can be filled arbitrarily
    putWord32le headerLength -- offset of pixel array in the file (bytes)

putBitmapInfoHeader (width, height) bpp = do
    putWord32le 40 -- size of BITMAPINFOHEADER
    putWord32le $ fromIntegral width
    putWord32le $ fromIntegral (-height) -- negative height -> first row is top of image
    putWord16le 1 -- # of color planes
    putWord16le bpp -- bits per pixel
    putWord32le 0 -- compression method (uncompressed RGB is 0)
    putWord32le 0 -- uncompressed size (optional if no compression is used)
    replicateM_ 2 $ putWord32le 0 -- dummy values for {x,y}-resolution
    putWord32le 0 -- use all 2^bpp colors
    putWord32le 0 -- # of "important colors", which is "generally ignored"

addBitmapHeader (w, h) image = do
    let headerLength = (2 + 4 + 4 + 4) + 40
    let imageLength = fromIntegral $ L.length image
    putBMPHeader headerLength imageLength
    putBitmapInfoHeader (w, h) 32
    putLazyByteString image

putRGB r g b = mapM_ putWord8 [b,g,r,0]

scaleFromTo (a, b) (a', b') x = ((x-a)/(b-a))*(b'-a')+a'

putPixel x y z = putRGB r g b where
    minimum = 0
    maximum = 2^8 - 1
    r = round $ if z < 0 then scaleFromTo (0, -1) (minimum, maximum) z else minimum
    g = round $ minimum
    b = round $ if z > 0 then scaleFromTo (0, 1) (minimum, maximum) z else minimum

visualizeFunction (xmin, ymin) (xmax, ymax) (w, h) (zmin, zmax) f = runPut image where
    xs = linspace w (xmin, xmax)
    ys = linspace h (ymin, ymax)
    f' = scaleFromTo (zmin, zmax) (-1, 1) . f
    image = V.forM_ ys (\y -> V.forM_ xs (\x -> putPixel x y (f' $ vector [x,y])))

visualizeFunctionBMP min max dims range = runPut . addBitmapHeader dims . visualizeFunction min max dims range

saveImage name img = withBinaryFile name WriteMode $ \h -> L.hPutStr h img

emitSample = saveImage "sample.bmp" $ visualizeFunctionBMP (-1,-1) (1,1) (2048, 2048) (-1,1) (tanh . V.sum)


text:
-- Bicubic interpolation algorithm from http://books.google.fi/books?id=1aAOdzK3FegC&pg=PA136&dq=bicubic+interpolation&hl=fi&sa=X&ei=japrUaWaH8Tcsga5hYHoCg&redir_esc=y#v=onepage&q=bicubic%20interpolation&f=false

module Gradient (gradientStepper) where

import Prelude hiding ((<>))
import Numeric.LinearAlgebra
import Types

bcuMatrix = (16><16)
  [
   1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
  -3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1, 0, 0, 0, 0,
   2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0,
   0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
   0, 0, 0, 0,-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1,
   0, 0, 0, 0, 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1,
  -3, 3, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0,
   9,-9, 9,-9, 6, 3,-3,-6, 6,-6,-3, 3, 4, 2, 1, 2,
  -6, 6,-6, 6,-4,-2, 2, 4,-3, 3, 3,-3,-2,-1,-1,-2,
   2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0,
  -6, 6,-6, 6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1,
   4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1::Double
  ]

--  0  1  2  3
--  4  5  6  7
--  8  9 10 11
-- 12 13 14 15
--
-- 01  (0,0) --> (0,1)   f, dx, dy, cross
-- 32    |    y=x1, x=x2
--       v
--     (1,0)
centeredDifferencing = (16><16) [
--    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15
      0,   0,   0,   0,   0,   1,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
      0,   0,   0,   0,   0,   0,   1,   0,   0,   0,   0,   0,   0,   0,   0,   0,
      0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   1,   0,   0,   0,   0,   0,
      0,   0,   0,   0,   0,   0,   0,   0,   0,   1,   0,   0,   0,   0,   0,   0,
      0,   0,   0,   0,-1/2,   0, 1/2,   0,   0,   0,   0,   0,   0,   0,   0,   0,
      0,   0,   0,   0,   0,-1/2,   0, 1/2,   0,   0,   0,   0,   0,   0,   0,   0,
      0,   0,   0,   0,   0,   0,   0,   0,   0,-1/2,   0, 1/2,   0,   0,   0,   0,
      0,   0,   0,   0,   0,   0,   0,   0,-1/2,   0, 1/2,   0,   0,   0,   0,   0,
      0,-1/2,   0,   0,   0,   0,   0,   0,   0, 1/2,   0,   0,   0,   0,   0,   0,
      0,   0,-1/2,   0,   0,   0,   0,   0,   0,   0, 1/2,   0,   0,   0,   0,   0,
      0,   0,   0,   0,   0,   0,-1/2,   0,   0,   0,   0,   0,   0,   0, 1/2,   0,
      0,   0,   0,   0,   0,-1/2,   0,   0,   0,   0,   0,   0,   0, 1/2,   0,   0,
      1/4, 0,-1/4,   0,   0,   0,   0,   0,-1/4,   0, 1/4,   0,   0,   0,   0,   0,
      0, 1/4,   0,-1/4,   0,   0,   0,   0,   0,-1/4,   0, 1/4,   0,   0,   0,   0,
      0,   0,   0,   0,   0, 1/4,   0,-1/4,   0,   0,   0,   0,   0,-1/4,   0, 1/4,
      0,   0,   0,   0, 1/4,   0,-1/4,   0,   0,   0,   0,   0,-1/4,   0, 1/4,   0::Double
  ]

combinedMatrix = bcuMatrix <> centeredDifferencing

-- bicubic interpolation
-- returns (value, d/dx, d/dy) at (y,x)@m
bicubicInterpolation :: Matrix Double -> (Double, Double) -> (Double, Double, Double)
bicubicInterpolation m (y, x) = (val, g1, g2)
  where
  (x_, dx) = properFraction x --u
  (y_, dy) = properFraction y --t
  ls  = tr $ reshape 4 $ combinedMatrix #> flatten sub
  us  = 4 |> [1.0,  dx,  dx*dx,  dx*dx*dx]
  ts  = 4 |> [1.0,  dy,  dy*dy,  dy*dy*dy]
  us2 = 4 |> [0.0, 1.0, 2.0*dx, 3.0*dx*dx]
  ts2 = 4 |> [0.0, 1.0, 2.0*dy, 3.0*dy*dy]
  val = ts  `dot` (ls #> us)  --f
  g1  = ts2 `dot` (ls #> us)  --d/dy
  g2  = ts  `dot` (ls #> us2) --d/dx
  sub = subMatrix (y_-1,x_-1) (4,4) m

gradientStepper:: Direction -> Double -> Matrix Double -> (Double, Double) -> [(Double, Double, Double)]
gradientStepper dir delta mat (y0, x0) = (x0, y0, z0) : ga x0 y0 z0 dy0 dx0
  where
  (z0, dy0, dx0) = bicubicInterpolation mat (y0, x0)
  ga x y z dy dx
    | x' < 2 || y' < 2 || x' > fromIntegral (cols mat - 3) || y' > fromIntegral (rows mat - 3) = []
    | z' > z = (x', y', z') : ga x' y' z' dy' dx'
    | otherwise = []
    where
    l
      |  dy == 00 && dx == 0.0 = 1
      | otherwise = 1.0 / sqrt (dy*dy+dx*dx)
    (dyn, dxn) = (dy*l, dx*l)
    directionOp = if dir == Descent then subtract else (+)
    (y',x') = (y `directionOp` delta*dyn, x `directionOp` delta*dxn)
    (z',dy',dx') = bicubicInterpolation mat (y', x')

max_stars_repo_name: angs/Hachure-map

max_issues_repo_name: angs/Hachure-map

max_forks_repo_name: angs/Hachure-map


text:
-- |

module DaySeven (examineFile) where
import Data.Vector.Unboxed (Vector, fromList, (!))
import qualified Data.Vector.Unboxed as V
import Statistics.Sample (mean)
import Statistics.Function (sort)

findOptimalFuelp2 :: Vector Int -> Int
findOptimalFuelp2 xs = minimum $ map (`calcFuel` xs) $ meanPositions xs
  where
    calcFuel = calcFuelBy sumOfNats

findOptimalFuelp1 :: Vector Int -> Int
findOptimalFuelp1 xs = (`calcFuel` xs) $ medianPos xs
  where
    calcFuel = calcFuelBy id

calcFuelBy :: (Int -> Int) -> Int -> Vector Int -> Int
calcFuelBy f x = V.foldl (\total pos -> total + f (abs (pos - x))) 0

-- | the sum of natural numbers
sumOfNats :: Int -> Int
sumOfNats n =  ((n*n) + n) `div` 2

meanPositions :: Vector Int -> [Int]
meanPositions vs = (map floor <> map ceiling) [mean $ V.map fromIntegral vs]

medianPos :: Vector Int -> Int
medianPos = floor . (\v -> sort v ! floor (fromIntegral (V.length v) / 2)) . V.map fromIntegral

examineFile :: String -> IO [Int]
examineFile xs = do
  positions <- fromList . (\x -> read $ "["++x++"]") <$> readFile xs
  return $ map (\f -> f positions) [findOptimalFuelp1,findOptimalFuelp2]


text:
-- Copyright (c) 1982-1999 Lennart Augustsson, Thomas Johnsson
-- See LICENSE for the full license.
--
module Hash(Hash(..), combineHash, emptyHash, hashToInt, hashToMax, Hashable(..)) where

--
-- Hash a value.  Hashing produces an Int of
-- unspecified range.
--
import Data.Array
import Data.Complex
import Data.Ratio

-- #define HASH(x) (case hash (x) of H h -> h)
hASH x = case hash (x) of H h -> h

newtype Hash = H Int deriving (Eq)

instance Show Hash where
    --showsType _ = showString "Hash"
    show _ = "<<Hash>>"

combineHash :: Hash -> Hash -> Hash
combineHash (H x) (H y) = H (x+y)

emptyHash :: Hash
emptyHash = H 0

class Hashable a where
    hash :: a -> Hash

instance Hashable Char where
    hash x = H $ fromEnum x

instance Hashable Int where
    hash x = H $ x

instance Hashable Integer where
    hash x = H $ fromInteger x

instance Hashable Float where
    hash x = H $ truncate x

instance Hashable Double where
    hash x = H $ truncate x

instance Hashable (IO a) where
    hash x = H 0

instance Hashable () where
    hash x = H 0

instance Hashable (a -> b) where
    hash x = H 0

instance (Hashable a) => Hashable (Maybe a) where
    hash Nothing = H 0
    hash (Just x) = hash x

instance (Hashable a, Hashable b) => Hashable (Either a b) where
    hash (Left x) = hash x
    hash (Right y) = hash y

instance Hashable a => Hashable [a] where
    hash l = H $ f l 0
        where f :: (Hashable a') => [a'] -> Int -> Int
              f [] r = r
              f (c:cs) r = f cs (3*r + hASH(c))

{-
#ifndef __HBCC__
instance Hashable [Char] where
    hash l = H $ f l 0
        where f :: String -> Int -> Int
              f [] r = r
              f (c:cs) r = f cs (3*r + fromEnum c)
#endif
-}

instance (Hashable a, Hashable b) => Hashable (a,b) where
    hash (a,b) = H $ hASH(a) + 3 * hASH(b)

instance (Hashable a, Hashable b, Hashable c) => Hashable (a,b,c) where
    hash (a,b,c) = H $ hASH(a) + 3 * hASH(b) + 5 * hASH(c)

instance (Hashable a, Hashable b, Hashable c, Hashable d) => Hashable (a,b,c,d) where
    hash (a,b,c,d) = H $ hASH(a) + 3 * hASH(b) + 5 * hASH(c) + 7 * hASH(d)

instance (Hashable a, Hashable b, Hashable c, Hashable d, Hashable e) => Hashable (a,b,c,d,e) where
    hash (a,b,c,d,e) = H $ hASH(a) + 3 * hASH(b) + 5 * hASH(c) + 7 * hASH(d) + 9 * hASH(e)

instance Hashable Bool where
    hash False = H 0
    hash True = H 1

instance (Integral a, Hashable a) => Hashable (Ratio a) where
    hash x = H $ hASH (denominator x) + hASH (numerator x)

instance (RealFloat a, Hashable a) => Hashable (Complex a) where
    hash (x :+ y) = H $ hASH(x) + hASH(y)

instance (Ix a) => Hashable (Array a b) where
    hash x = H $ 0 -- !!!

hashToInt :: Int -> Hash -> Int
hashToInt maxhash x =
    case x of
    H h ->
        if h < 0 then
            if -h < 0 then 0
            else (-h) `rem` maxhash
        else h `rem` maxhash

hashToMax maxhash x =
    case hash x of
    H h ->
        if h < 0 then
            if -h < 0 then 0
            else (-h) `rem` maxhash
        else h `rem` maxhash

max_stars_repo_name: GaloisInc/BESSPIN-BSC

max_issues_repo_name: GaloisInc/BESSPIN-BSC

max_forks_repo_name: GaloisInc/BESSPIN-BSC


text:
import Data.Complex (cis, phase)

meanAngle
  :: RealFloat c
  => [c] -> c
meanAngle = (/ pi) . (* 180) . phase . sum . map (cis . (/ 180) . (* pi))

main :: IO ()
main =
  mapM_
    (\angles ->
        putStrLn $
        "The mean angle of " ++
        show angles ++ " is: " ++ show (meanAngle angles) ++ " degrees")
    [[350, 10], [90, 180, 270, 360], [10, 20, 30]]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE NamedFieldPuns #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}

import qualified Data.Map as Map
import Data.Map(Map)
import qualified Data.Set as Set
import Control.Arrow(first, second, (&&&))
import Data.List
import Data.Maybe
import Debug.Trace
import Data.Ord
import qualified Text.Printf as Printf
import qualified Data.Array as Array
import Data.Array(Array, (!), range, Ix)
import qualified Data.Matrix
import qualified System.Random
import Control.Monad
import Control.Monad.Identity
import System.IO(stdout, hSetBuffering, BufferMode(..))
import Control.Parallel.Strategies
import GHC.Generics
import qualified Data.Eigen.Matrix
import Foreign.C.Types
import qualified Numeric.LinearAlgebra.HMatrix
import qualified Numeric.LinearAlgebra.Data
import Control.Arrow
import GHC.Stack

import Matrix

data RawMaterial =
  BuriedIronRaw
  | BuriedCopperRaw
  | BuriedStoneRaw
--  | PetroleumGasRaw
--  | LightOilRaw
--  | HeavyOilRaw
  | CrudeOilRaw
  | BuriedCoalRaw
  | PollutionRaw
  deriving (Eq, Ord, Enum, Bounded, Show, Ix, Generic)

instance NFData RawMaterial

type RawMaterialPressure = Map Product Rat

data Product =
  IronPlate
  | CopperPlate
  | SteelPlate
  | IronOre
  | CopperOre
  | CopperCable
  | GearWheel
  | Plastic
  | Pipe
  | ElectronicCircuit
  | AdvancedCircuit
  | EngineUnit
  | ProcessingUnit
  | ElectricMiningDrill
  | PiercingRoundMagazine
  | FirearmRoundMagazine
  | Grenade
  | GunTurret
  | SciencePack1
  | SciencePack2
  | SciencePack3
  | SciencePackMilitary
  | SciencePackProduction
  | SciencePackHighTech
  | SpeedModule
  | SpeedModule2
  | SpeedModule3
  | ProductivityModule
  | ProductivityModule2
  | ProductivityModule3
  | EfficiencyModule
  | EfficiencyModule2
  | EfficiencyModule3
  | Sulfur
  | SulfuricAcid
  | CoalLiquefaction
  | ResearchCoalLiquefaction
  | ResearchLaserTurretDamage5
  | ResearchRocketSilo
  | ResearchNuclearPower
  | ResearchEndgame
  | Inserter
  | TransportBelt
  | AssemblingMachine1
  | AssemblingMachine2
  | AssemblingMachine3
  | SteelFurnace
  | ElectricFurnace
  | ElectricEngineUnit
  | Lubricant
  | StoneBrick
  | Stone
  | LaserTurret
  | Battery
  | ChemicalPlant
  | OilRefinery
  | LabBuilding
  | BoilerBuilding
  | SteamEngineBuilding
  | SolarFacilityBuilding
  | NuclearFacilityBuilding
  | Roboport
  | Substation
  | Accumulator
  | SolarPanel
  | StoneFurnace
  | SteamTurbine
  | NuclearReactor
  | HeatExchanger
  | HeatPipe
  | Concrete
  | Beacon
  | RocketPart
  | RocketSiloBuilding
  | ControlModule
  | LightweightStructure
  | RocketFuel
  | SciencePackSpace
  | Satellite
  | Radar

  | ElectricalEnergy -- in J
  | ChemicalEnergy -- in J
  | SolidFuel
  | Steam
  | Coal
  | BuriedIron
  | BuriedCopper
  | PetroleumGas
  | LightOil
  | HeavyOil
  | CrudeOil
  | BuriedCoal
  | BuriedStone
  | Pollution
  deriving (Eq, Ord, Show, Enum, Bounded, Ix, Generic)

instance NFData Product where

newtype Time = Time { unTime :: Rat } deriving (Eq, Ord, Show, NFData)

instance Linear Time where
  zero = Time 0
  add (Time x) (Time y) = Time (add x y)
  minus (Time x) = Time (minus x)

instance VectorSpace Time where
  type Scalar Time = Rat
  scale x (Time t) = Time (x * t)

data VenueKind =
  AssemblyVenueKind
  | SmelterVenueKind
  | ChemicalVenueKind
  | MinerVenueKind
  | LabVenueKind
  | BoilerVenueKind
  | SteamEngineVenueKind
  | GreenPowerVenueKind
  | RefineryVenueKind
  | RocketSiloVenueKind
  | NoVenueVenueKind
  deriving (Show, Ord, Eq)
venuesByKind :: VenueKind -> [Venue]
venuesByKind AssemblyVenueKind = [Assembly2, Assembly3]
venuesByKind SmelterVenueKind = [SmelterBurner, SmelterElectric]
venuesByKind ChemicalVenueKind = [Chemical]
venuesByKind MinerVenueKind = [Miner]
venuesByKind LabVenueKind = [Lab]
venuesByKind BoilerVenueKind = [Boiler]
venuesByKind SteamEngineVenueKind = [SteamEngine]
venuesByKind GreenPowerVenueKind = [SolarFacility, NuclearFacility]
venuesByKind RefineryVenueKind = [Refinery]
venuesByKind NoVenueVenueKind = [NoVenue]
venuesByKind RocketSiloVenueKind = [RocketSilo]

data Venue =
  Assembly2
  | Assembly3
  | SmelterBurner
  | SmelterElectric
  | Chemical
  | Miner
  | Lab
  | Boiler
  | SteamEngine
  | SolarFacility
  | NuclearFacility
  | Refinery
  | RocketSilo
  | NoVenue
  deriving (Show, Eq, Ord, Generic)

data ModuleConfig = ModuleConfig
  {
    configSpeedBonus :: Rat,
    configProductivityBonus :: Rat,
    configEnergyBonus :: Rat,
    configPollutionBonus :: Rat
  } deriving (Eq, Ord, Show, Generic)

instance Linear ModuleConfig where
  zero = mempty
  add = mappend
  minus = scale (-1)

instance VectorSpace ModuleConfig where
  type Scalar ModuleConfig = Rat
  scale x (ModuleConfig a b c d) =
    ModuleConfig
     (x * a) (x * b) (x * c) (x * d)

-- venue, modules and modules in beacons
type PreConfig = (Venue, [Product], Maybe Product)

data Config = Config
  {
    configVenue :: Venue,
    configModules :: ModuleConfig,
    configModuleMaterials :: Map Product Rat,
    configConstantExtraPower :: Rat
  } deriving Generic

instance NFData ModuleConfig
instance NFData Config

speedMultiplier c = (1 + configSpeedBonus x) where x = configModules c
productivityMultiplier c = 1 + configProductivityBonus x where x = configModules c
energyMultiplier c = max 0.2 (1 + configEnergyBonus x) where x = configModules c
pollutionMultiplier c = (1 + configPollutionBonus x) where x = configModules c

instance Monoid ModuleConfig where
  mempty = ModuleConfig 0 0 0 0
  a `mappend` b =
    ModuleConfig
      {
        configSpeedBonus = (configSpeedBonus a + configSpeedBonus b),
        configProductivityBonus = (configProductivityBonus a + configProductivityBonus b),
        configEnergyBonus = (configEnergyBonus a + configEnergyBonus b),
        configPollutionBonus = (configPollutionBonus a + configPollutionBonus b)
      }

data Usability =
  Unusable | Usable

moduleToConfig SpeedModule = ModuleConfig 0.2 0 0.5 0
moduleToConfig SpeedModule2 = ModuleConfig 0.3 0 0.6 0
moduleToConfig SpeedModule3 = ModuleConfig 0.5 0 0.7 0
moduleToConfig EfficiencyModule = ModuleConfig 0 0 (negate 0.3) 0
moduleToConfig EfficiencyModule2 = ModuleConfig 0 0 (negate 0.4) 0
moduleToConfig EfficiencyModule3 = ModuleConfig 0 0 (negate 0.5) 0
moduleToConfig ProductivityModule = ModuleConfig (negate 0.15) 0.04 0.4 0
moduleToConfig ProductivityModule2 = ModuleConfig (negate 0.15) 0.06 0.6 0
moduleToConfig ProductivityModule3 = ModuleConfig (negate 0.15) 0.10 0.8 0
moduleToConfig p = error $ "not a module: " ++ show p

allModules :: Usability -> [[[Product]]]
allModules usability =
  [[ ([]) ]] ++
  (map (\ms -> map (\m -> [m]) ms) $
  (map return
  [ SpeedModule
  , EfficiencyModule
--  , EfficiencyModule2
--  , EfficiencyModule3
  , SpeedModule2
  , SpeedModule3
  ]
  ++ [case usability of
   Unusable ->
     [ ProductivityModule
     , ProductivityModule2
     , ProductivityModule3
     ]
   Usable -> []]))

choose' :: Int -> [[a]] -> [[a]]
choose' k l
  | k < 0 = []
  | k == 0 = [[]]
  | otherwise = case l of
    [] -> []
    (xc : xs) -> choose' k xs ++ (do
      i <- [1..k]
      (x :: a) <- xc
      rest <- choose' (k - i) xs
      return $ replicate i x ++ rest)

choose k l
  | k < 0 = []
  | k == 0 = [[]]
  | otherwise = case l of
     [] -> []
     (x : xs) -> map (x:) (choose (k-1) l) ++ choose k xs

moduleSlots :: Venue -> Int
moduleSlots venue = case venue of
        Assembly2 -> 2
        Assembly3 -> 4
        SmelterElectric -> 2
        SmelterBurner -> 0
        Chemical -> 3
        Miner -> 3
        Lab -> 2
        Boiler -> 0
        SteamEngine -> 0
        SolarFacility -> 0
        NuclearFacility -> 0
        Refinery -> 3
        NoVenue -> 0
        RocketSilo -> 4

initial_module_config gc Miner = mempty { configProductivityBonus = gc_miner_productivity_bonus gc }
initial_module_config _ _ = mempty

-- we have [beaconizationFactor] factories per beacon
-- and [beaconizationFactor] beacons per factory
beaconizationFactor :: RecipeName -> Integer
beaconizationFactor (ProductRecipe ElectronicCircuit) = 2
beaconizationFactor (ProductRecipe CopperCable) = 7
beaconizationFactor (ProductRecipe SciencePackProduction) = 1
beaconizationFactor (ProductRecipe SciencePackMilitary) = 2
beaconizationFactor (ProductRecipe IronOre) = 6
beaconizationFactor _ = 4

linksPerBeacon :: RecipeName -> Integer
linksPerBeacon = beaconizationFactor
linksPerFactory :: RecipeName -> Integer
linksPerFactory = beaconizationFactor

beaconsPerFactory x = fromInteger (linksPerFactory x) / fromInteger (linksPerBeacon x)

mkConfig gc recipe (venue, modules, beacon) =
  Config
    venue
    (mconcat $
      (initial_module_config gc venue)
      : map moduleToConfig modules
      ++ (case beacon of
        Just m ->
          [ scale (fromIntegral (linksPerFactory recipe)) (moduleToConfig m) ]
        Nothing ->
          []
        )
      )
    (Map.fromListWith add
    (map (\x -> (x, 1)) modules ++ (
      case beacon of
       Nothing -> []
       Just m ->
         let ratio = beaconsPerFactory recipe in
         [(m, 2 * ratio), (Beacon, 1 * ratio)]
     )
    ))
    (case beacon of
      Nothing -> 0
      Just _ -> 480e3 * (beaconsPerFactory recipe)
    )

-- in Watt
data Power =
  ElectricalPower { unPower :: Rat }
  | ChemicalPower Rat

baseBoilerPower = 3.6e6

basePower :: Venue -> Power
basePower Assembly2 = ElectricalPower 150e3
basePower Assembly3 = ElectricalPower 210e3
basePower Miner = ElectricalPower 90e3
basePower SmelterElectric = ElectricalPower 180e3
basePower SmelterBurner = ChemicalPower 180e3
basePower Chemical = ElectricalPower 210e3
basePower Lab = ElectricalPower 60e3
basePower Boiler = ChemicalPower baseBoilerPower
basePower SteamEngine = ElectricalPower 0
basePower SolarFacility = ElectricalPower 0
basePower NuclearFacility = ElectricalPower 0
basePower Refinery = ElectricalPower 420e3
basePower NoVenue = ElectricalPower 0
basePower RocketSilo = ElectricalPower 4e6

basePollution Assembly2 = 2.4
basePollution Assembly3 = 1.8
basePollution Miner = 9
basePollution SmelterElectric = 0.9
basePollution SmelterBurner = 3.6
basePollution Chemical = 1.8
basePollution Lab = 0
basePollution Boiler = 27.6923
basePollution Refinery = 3.6
basePollution SolarFacility = 0
basePollution NuclearFacility = 0
basePollution SteamEngine = 0
basePollution NoVenue = 0
basePollution RocketSilo = 0

baseSpeed _ SolarFacility = (42e3 * 176)
baseSpeed _ NuclearFacility = 160e6
baseSpeed _ Assembly2 = 0.75
baseSpeed _ Assembly3 = 1.25
baseSpeed _ Miner = 1 -- factored in into the recipe
baseSpeed _ SmelterElectric = 2
baseSpeed _ SmelterBurner = 2
baseSpeed _ Chemical = 1.25
baseSpeed config Lab = gc_lab_speed_multiplier config
baseSpeed _ Boiler = 1 -- this is factored into the recipe
baseSpeed _ SteamEngine = 1 -- this is factored into the recipe
baseSpeed _ Refinery = 1
baseSpeed _ NoVenue = 1 -- this is meaningless
baseSpeed _ RocketSilo = 1

scaleTime s (Time t) = Time (s * t)

coalToEnergy coal = coal * 8e6 / 2

data Recipe = Recipe
  {
    recipeName :: RecipeName,
    recipeProducts :: [(Product, Rat)],
    recipeMaterials :: [(Product, Rat)],
    recipeVenueKind :: VenueKind,
    recipeTime :: Time
  } deriving (Eq, Ord, Show)

data RecipeName =
  ProductRecipe Product
  | LiquefactionRecipe
  | AdvancedOilProcessing
  | BoilerRecipe
  | UseAsFuelRecipe Product
  deriving (Eq, Ord, Generic)

instance Show RecipeName where
  show (ProductRecipe product) = show product
  show (LiquefactionRecipe) = "Liquefaction"
  show (AdvancedOilProcessing) = "AdvancedOilProcessing"
  show BoilerRecipe = "BoilerRecipe"
  show (UseAsFuelRecipe product) = "UseAsFuel" ++ show product

energy_per_steam = 30000

data ResearchCost = ResearchCost {
  researchCostTime :: Time, -- time to process 1 stack of bottles
  researchCostBottles :: [Product],
  researchCostAmount :: Rat
}

marathon_research_cost_adjustment = 4.0

lab_speed_researches :: [(ResearchCost, Rat)]
lab_speed_researches =
  let r = SciencePack1 in
  let g = SciencePack2 in
  let b = SciencePack3 in
  let p = SciencePackProduction in
  let y = SciencePackHighTech in
  let bot l amount =  ResearchCost (Time 30) l (amount * marathon_research_cost_adjustment) in
  (bot [r, g] 100, 0.20)
  : (bot [r, g] 200, 0.30)
  : (bot [r, g, b] 250, 0.40)
  : (bot [r, g, b] 500, 0.50)
  : (bot [r, g, b, p] 500, 0.50)
  : (bot [r, g, b, p, y] 500, 0.60)
  : error "further lab researches unknown"

mining_productivity_researches :: [(ResearchCost, Rat)]
mining_productivity_researches =
  let r = SciencePack1 in
  let g = SciencePack2 in
  let b = SciencePack3 in
  let p = SciencePackProduction in
  let y = SciencePackHighTech in
  zipWith
    (\colors amounts -> (ResearchCost (Time 60) colors (amounts * 4), 0.02))
    (replicate 3 [r, g] ++ replicate 4 [r, g, b] ++ replicate 4 [r, g, b, p] ++ replicate 4 [r, g, b, p, y])
    (map fromInteger [100, 200..]) ++ error "further mining productivity unknown"

computeResearchCapital :: Int -> [(ResearchCost, Rat)] -> Map Product Rat
computeResearchCapital n =
  mconcat' . map (\r -> scale (researchCostAmount r) (Map.fromListWith (+) $ map (\p -> (p, 1)) (researchCostBottles r))) . map fst . take n

data SmeltingMode =
  BurnerSmelting
  | ElectricSmelting

data Liquefaction_mode =
  Liquefaction_disabled
  | Liquefaction_to_gas
  | Liquefaction_to_burn
   deriving (Show, Generic)

class Enumerate a where
  allOfThem :: [a]

instance Enumerate Liquefaction_mode where
  allOfThem =
   [ Liquefaction_disabled
   , Liquefaction_to_gas
   , Liquefaction_to_burn
   ]

-- the details of GameConfig sufficient to determine what recipes are available
data GameConfigQualitative = GameConfigQualitative {
  qgc_liquefaction :: Liquefaction_mode
  }

to_qualitative (GameConfig { gc_liquefaction = liquefaction }) = GameConfigQualitative { qgc_liquefaction = liquefaction }

data GameConfig = GameConfig {
  gc_lab_researches_done :: Int,
  gc_mining_researches_done :: Int,
  gc_liquefaction :: Liquefaction_mode,
  gc_recipe_configs :: RecipeName -> PreConfig
  }
  deriving Generic

data Change =
  ProductChange RecipeName PreConfig
  | Other String String
  deriving Generic

instance NFData Change
instance NFData GameConfig
instance NFData Liquefaction_mode
instance NFData RecipeName

gc_alternatives :: GameConfig -> [(Change, GameConfig)]
gc_alternatives gc =
  [(Other "Liquefaction Mode" (show mode), gc { gc_liquefaction = mode }) | mode <- allOfThem] ++
  [(Other "Lab speed" "+1", gc { gc_lab_researches_done = gc_lab_researches_done gc + 1 }) ] ++
  [(Other "Mining productivity" "+1", gc { gc_mining_researches_done = gc_mining_researches_done gc + 1 }) ] ++
  [(ProductChange recipeName config, gc { gc_recipe_configs = f' })
  | (recipe, _) <- recipes
  , recipeName <- [ recipeName recipe ]
  , venue <- venuesByKind (recipeVenueKind recipe)
  , let availableModules = allModules (usability recipeName)
  , modules <- choose' (moduleSlots venue) availableModules
  , modules <- [ concat modules ]
  , beacon <- if length modules == moduleSlots venue && moduleSlots venue > 0 then [ Nothing, Just SpeedModule2, Just SpeedModule3 ] else [ Nothing ]
  , let config = (venue, modules, beacon)
  , let f' = (let f = gc_recipe_configs gc in (\r -> if r == recipeName then config else f r))
  ]

gc_lab_speed_multiplier :: GameConfig -> Rat
gc_lab_speed_multiplier t = 1 +
  (mconcat' $ map snd $ take (gc_lab_researches_done t) lab_speed_researches)

gc_miner_productivity_bonus :: GameConfig -> Rat
gc_miner_productivity_bonus t =
  (mconcat' $ map snd $ take (gc_mining_researches_done t) mining_productivity_researches)

gc_configs :: GameConfig -> (Recipe -> Config)
gc_configs gc =
 \recipe ->
 let pre = gc_recipe_configs gc (recipeName recipe) in
 mkConfig gc (recipeName recipe) pre

-- recipes and whether or not they should be enabled
recipes :: [(Recipe, (GameConfigQualitative -> Bool))]
recipes =
  let assembly = AssemblyVenueKind in
  let smelter = SmelterVenueKind in
  let
   miner what ingredients speed =
     [ (Recipe (ProductRecipe what) [(what, 1)] ingredients MinerVenueKind (Time $ recip speed), const True) ]
  in
  concat [
    r GearWheel 1 [(IronPlate, 4)] assembly (Time 0.5),
    r IronPlate 1 [(IronOre, 1)] smelter (Time 3.5),
    r CopperPlate 1 [(CopperOre, 1)] smelter (Time 3.5),
    r SteelPlate 1 [(IronPlate, 10)] smelter (Time 35),
    miner IronOre [(BuriedIron, 1)] 0.525,
    miner CopperOre [(BuriedCopper, 1)] 0.525,
    miner Coal [(BuriedCoal, 1)] 0.525,
    r Beacon 1[(AdvancedCircuit, 20), (CopperCable, 10), (ElectronicCircuit, 20), (SteelPlate, 10)] assembly (Time 15),
    r Plastic 2[(PetroleumGas, 20), (Coal, 1)] ChemicalVenueKind (Time 1),
    r ElectronicCircuit 1 [(CopperCable, 10), (IronPlate, 2)] assembly (Time 0.5),
    r AdvancedCircuit 1 [(Plastic, 4), (CopperCable, 8), (ElectronicCircuit, 2)] assembly (Time 6),
    r CopperCable 2 [(CopperPlate, 1)] assembly (Time 0.5),
    r Pipe 1 [(IronPlate, 2)] assembly (Time 0.5),
    r EngineUnit 1 [(GearWheel, 1), (Pipe, 2), (SteelPlate, 1)] assembly (Time 10),
    r ElectricMiningDrill 1 [(GearWheel, 10), (IronPlate, 20), (ElectronicCircuit, 5)] assembly (Time 2),
    r SciencePack3 1 [(AdvancedCircuit, 1), (ElectricMiningDrill, 1), (EngineUnit, 1)] assembly (Time 12),
    r SciencePack1 1 [(CopperPlate, 1), (GearWheel, 1)] assembly (Time 5),
    r SciencePack2 1 [(Inserter, 1), (TransportBelt, 1)] assembly (Time 6),
    r SciencePackProduction 2 [(ElectricEngineUnit, 1), (ElectricFurnace, 1)] assembly (Time 14),
    r PiercingRoundMagazine 1 [(FirearmRoundMagazine, 1), (SteelPlate, 1), (CopperPlate, 5)] assembly (Time 3),
    r FirearmRoundMagazine 1 [(IronPlate, 4)] assembly (Time 1),
    r Grenade 1 [(IronPlate, 5), (Coal, 10)] assembly (Time 8),
    r GunTurret 1[(GearWheel, 10), (CopperPlate, 10), (IronPlate, 20)] assembly (Time 8),
    r SciencePackMilitary 2 [(PiercingRoundMagazine, 1), (Grenade, 1), (GunTurret, 1)] assembly (Time 10),
    r SpeedModule 1 [(AdvancedCircuit, 5), (ElectronicCircuit, 5)] assembly (Time 15),
    r EfficiencyModule 1 [(AdvancedCircuit, 5), (ElectronicCircuit, 5)] assembly (Time 15),
    r ProductivityModule 1 [(AdvancedCircuit, 5), (ElectronicCircuit, 5)] assembly (Time 15),
    r EfficiencyModule2 1 [(AdvancedCircuit, 5), (EfficiencyModule, 4), (ProcessingUnit, 5)] assembly (Time 30),
    r SpeedModule2 1 [(AdvancedCircuit, 5), (SpeedModule, 4), (ProcessingUnit, 5)] assembly (Time 30),
    r ProductivityModule2 1 [(AdvancedCircuit, 5), (ProductivityModule, 4), (ProcessingUnit, 5)] assembly (Time 30),
    r ProductivityModule3 1 [(AdvancedCircuit, 5), (ProductivityModule2, 5), (ProcessingUnit, 5)] assembly (Time 60),
    r SpeedModule3 1 [(AdvancedCircuit, 5), (SpeedModule2, 5), (ProcessingUnit, 5)] assembly (Time 60),
    r EfficiencyModule3 1 [(AdvancedCircuit, 5), (EfficiencyModule2, 5), (ProcessingUnit, 5)] assembly (Time 60),
    r ProcessingUnit 1 [(AdvancedCircuit, 2), (ElectronicCircuit, 20), (SulfuricAcid, 10)] assembly (Time 10),
    r SulfuricAcid 50 [(IronPlate, 1), (Sulfur, 5)] ChemicalVenueKind (Time 1),
    r Sulfur 2 [(PetroleumGas, 30)] ChemicalVenueKind (Time 1),
    r ResearchCoalLiquefaction (1/800) [(SciencePack1, 1), (SciencePack2, 1), (SciencePack3, 1), (SciencePackProduction, 1)] LabVenueKind (Time 30),
    r ResearchNuclearPower (1/4000) [(SciencePack1, 1), (SciencePack2, 1), (SciencePack3, 1)] LabVenueKind (Time 30),
    r ResearchLaserTurretDamage5 (1/800)
      [(SciencePack1, 1), (SciencePack2, 1), (SciencePack3, 1), (SciencePackProduction, 1), (SciencePackHighTech, 1)] LabVenueKind (Time 60),
    r ResearchRocketSilo (1/4000)
      [ (SciencePack1, 1)
      , (SciencePack2, 1)
      , (SciencePack3, 1)
      , (SciencePackProduction, 1)
      , (SciencePackHighTech, 1)
      , (SciencePackMilitary, 1)] LabVenueKind (Time 60),
    r ResearchEndgame (1/10000)
      [ (SciencePack1, 1)
      , (SciencePack2, 1)
      , (SciencePack3, 1)
      , (SciencePackProduction, 1)
      , (SciencePackHighTech, 1)
      , (SciencePackMilitary, 1)
      , (SciencePackSpace, 1)
      ] LabVenueKind (Time 60),
    r Inserter 1 [(ElectronicCircuit, 1), (IronPlate, 1), (GearWheel, 1)] assembly (Time 0.5),
    r TransportBelt 2 [(GearWheel, 1), (IronPlate, 1)] assembly (Time 0.5),
    r AssemblingMachine1 1 [(GearWheel, 5), (IronPlate, 9), (ElectronicCircuit, 3)] assembly (Time 0.5),
    r AssemblingMachine2 1 [(AssemblingMachine1, 1), (ElectronicCircuit, 5), (GearWheel, 10), (IronPlate, 20)] assembly (Time 0.5),
    r AssemblingMachine3 1 [(AssemblingMachine2, 2), (SpeedModule, 4)] assembly (Time 0.5),
    r SteelFurnace 1 [(SteelPlate, 6), (StoneBrick, 10)] assembly (Time 3),
    r ElectricFurnace 1 [(AdvancedCircuit, 5), (SteelPlate, 10), (StoneBrick, 10)] assembly (Time 5),
    r ElectricEngineUnit 1 [(ElectronicCircuit, 2), (EngineUnit, 1), (Lubricant, 15)] assembly (Time 10),
    r StoneBrick 1 [(Stone, 2)] smelter (Time 3.5),
    miner Stone [(BuriedStone, 1)] 0.65,
    r Lubricant 10 [(HeavyOil, 10)] ChemicalVenueKind (Time 1),
    r LaserTurret 1 [(Battery, 12), (ElectronicCircuit, 20), (SteelPlate, 20)] assembly (Time 20),
    r Battery 1[(CopperPlate, 1), (IronPlate, 1), (SulfuricAcid, 40)] ChemicalVenueKind (Time 5),
    r SciencePackHighTech 2 [(Battery, 1), (CopperCable, 30), (ProcessingUnit, 3), (SpeedModule, 1)] assembly (Time 14),
    r ChemicalPlant 1 [(ElectronicCircuit, 5), (GearWheel, 5), (Pipe, 5), (SteelPlate, 5)] assembly (Time 5),
    r OilRefinery 1 [(ElectronicCircuit, 10), (GearWheel, 10), (Pipe, 10), (SteelPlate, 15), (StoneBrick, 10)] assembly (Time 8),
    r LabBuilding 1 [(ElectronicCircuit, 10), (GearWheel, 10), (TransportBelt, 4)] assembly (Time 2),
    r BoilerBuilding 1 [(Pipe, 4), (StoneFurnace, 1)] assembly (Time 0.5),
    r StoneFurnace 1 [(Stone, 5)] assembly (Time 0.5),
    r SteamEngineBuilding 1 [(GearWheel, 10), (IronPlate, 50), (Pipe, 5)] assembly (Time 0.5),

    [ (Recipe BoilerRecipe [(Steam, (baseBoilerPower * 0.5) / energy_per_steam)] [] BoilerVenueKind (Time 1), const True) ],
--    r ElectricalEnergy 1 [(Steam, 1/energy_per_steam)] SteamEngineVenueKind (Time (1/900e3)),
    r ElectricalEnergy 1 [] GreenPowerVenueKind (Time 1),
    r SolarFacilityBuilding 1 [(SolarPanel, 176), (Accumulator, 166), (Substation, 10), (Roboport, 1)] NoVenueVenueKind (Time 1),
    r Roboport 1 [(AdvancedCircuit, 45), (GearWheel, 45), (SteelPlate, 45)] assembly (Time 5),
    r Substation 1 [(AdvancedCircuit, 5), (CopperPlate, 5), (SteelPlate, 10)] assembly (Time 0.5),
    r Accumulator 1 [(Battery, 5), (IronPlate, 2)] assembly (Time 10),
    r SolarPanel 1 [(CopperPlate, 5), (ElectronicCircuit, 15), (SteelPlate, 5)] assembly (Time 10),

    r PetroleumGas 2 [(LightOil, 3)] ChemicalVenueKind (Time 5),

    r SolidFuel 1 [(LightOil, 10)] ChemicalVenueKind (Time 3),
    r LightOil 3 [(HeavyOil, 4)] ChemicalVenueKind (Time 5),
    [ (Recipe AdvancedOilProcessing [(HeavyOil, 10), (LightOil, 45), (PetroleumGas, 55)] [(CrudeOil, 100)] RefineryVenueKind (Time 5), const True) ],
    r Concrete 10 [(StoneBrick, 5), (IronOre, 1)] assembly (Time 10),
    r NuclearFacilityBuilding 1 [(NuclearReactor, 1), (HeatExchanger, 16), (SteamTurbine, 25), (HeatPipe, 60)] NoVenueVenueKind (Time 1),
    r SteamTurbine 1 [(GearWheel, 50), (CopperPlate, 50), (Pipe, 20)] assembly (Time 3),
    r NuclearReactor 1 [(Concrete, 500), (SteelPlate, 500), (AdvancedCircuit, 500), (CopperPlate, 500)] assembly (Time 8),
    r HeatExchanger 1 [(SteelPlate, 10), (CopperPlate, 100), (Pipe, 10)] assembly (Time 3),
    r HeatPipe 1 [(SteelPlate, 10), (CopperPlate, 20)] assembly (Time 1),

    r RocketSiloBuilding 1 [(Concrete, 1000), (ElectricEngineUnit, 200), (Pipe, 100), (ProcessingUnit, 200), (SteelPlate, 1000)] assembly (Time 30),
    r RocketPart 1 [(ControlModule, 10), (LightweightStructure, 10), (RocketFuel, 10)] RocketSiloVenueKind (Time 3),
    r ControlModule 1 [(ProcessingUnit, 1), (SpeedModule, 1)] assembly (Time 30),
    r LightweightStructure 1 [(CopperPlate, 10), (Plastic, 10), (SteelPlate, 10)] assembly (Time 30),
    r RocketFuel 1 [(SolidFuel, 10)] ChemicalVenueKind (Time 30),
    r SciencePackSpace 1000 [(RocketPart, 100), (Satellite, 1)] NoVenueVenueKind (Time 1),
    r Satellite 1 [(Accumulator, 100), (LightweightStructure, 100), (ProcessingUnit, 100), (Radar, 5), (RocketFuel, 50), (SolarPanel, 100)] assembly (Time 5),
    r Radar 1 [(ElectronicCircuit, 5), (GearWheel, 5), (IronPlate, 10)] assembly (Time 0.5),
    (
    let
     burn_coal =
      Recipe (UseAsFuelRecipe Coal) [(ChemicalEnergy, 8e6)] [(Coal, 1)] NoVenueVenueKind (Time 1) -- time is meaningless here
     liquefaction =
      Recipe LiquefactionRecipe [(HeavyOil, 35), (LightOil, 15), (PetroleumGas, 20)] [(Coal, 10), (HeavyOil, 25), (Steam, 50)] RefineryVenueKind (Time 5)
     burn_bricks =
      Recipe (UseAsFuelRecipe SolidFuel) [(Steam, (25e6 * 0.5) / energy_per_steam)] [(SolidFuel, 1)] BoilerVenueKind (Time 1) -- incorrect time, but nothing cares
    in
    [ (burn_coal, (\gc -> case qgc_liquefaction gc of { Liquefaction_to_burn -> False; _ -> True }))
    , (burn_bricks, (\gc -> case qgc_liquefaction gc of { Liquefaction_to_burn -> True; _ -> False }))
    , (liquefaction, (\gc -> case qgc_liquefaction gc of { Liquefaction_disabled -> False; _ -> True }))
    ])
  ] where
  r product quantity ingredients venues time = [ (Recipe (ProductRecipe product) [(product, quantity)] ingredients venues time, const True) ]

mconcat' x = foldr add zero x

functionToMatrix :: (Num v, Ix' b, Ix' a) => (a -> [(b, v)]) -> Matrix a b v
functionToMatrix f =
  Matrix (Array.array fullRange [ ((a, b), maybe 0 id (Map.lookup b bs)) | a <- range fullRange, let bs = Map.fromListWith (+) (f a), b <- range fullRange])

recipesByName = Map.fromListWith (error "multiple recipes with the same name") (map (\(recipe, _enabled) -> (recipeName recipe, recipe)) (recipes))

enabledRecipesByName gc =
  Map.fromListWith (error "multiple recipes with the same name") (concatMap (\(recipe, enabled) -> if not (enabled gc) then [] else [(recipeName recipe, recipe)]) (recipes))

recipesToMatrix :: GameConfig -> Map RecipeName (Map Product Rat)
recipesToMatrix gc =
  fmap (\(Recipe recipeName production consumption _venueKind (Time baseTime)) ->
         (
         let preConfig = gc_recipe_configs gc recipeName in
         let config = mkConfig gc recipeName preConfig in
         let venue = configVenue config in
           let
             time = (baseTime / (speedMultiplier config * baseSpeed gc venue))
             extra_energy = time * configConstantExtraPower config
             energy_and_pollution =
               let multiplier = time * energyMultiplier config in
               let pollution = basePollution venue in
               case basePower venue of
                 ElectricalPower basePower ->
                   [(ElectricalEnergy, (basePower * multiplier) + extra_energy), (Pollution, pollution * multiplier)]
                 ChemicalPower basePower ->
                   [(ChemicalEnergy, basePower * multiplier + extra_energy), (Pollution, pollution * multiplier)]
           in
             Map.fromListWith add (fmap (second negate) (consumption ++ energy_and_pollution) ++ map (second ((* productivityMultiplier config))) production)
         )) (enabledRecipesByName (to_qualitative gc))

matrix_add (Matrix x) (Matrix y) = Matrix (Array.array fullRange [(i, x ! i + y ! i) | i <- range fullRange])
matrix_negate (Matrix x) = Matrix (Array.array fullRange [(i, - x ! i) | i <- range fullRange])

matrix_subtract a b = matrix_add a (matrix_negate b)

matrix_identity :: (Ix' a, Num v) => Matrix a a v
matrix_identity = Matrix (Array.array fullRange [((a,b), if a == b then 1 else 0) | (a,b) <- range fullRange])

to_data_matrix :: forall a v . (Ix' a) => Matrix a a v -> Data.Matrix.Matrix v
to_data_matrix (Matrix m) =
  Data.Matrix.matrix
  (Array.rangeSize (fullRange :: (a, a)))
  (Array.rangeSize (fullRange :: (a, a)))
  (\(i, j) ->
     m ! (toEnum (i - 1), toEnum (j - 1)))

of_data_matrix :: forall a v . Ix' a => Data.Matrix.Matrix v -> Matrix a a v
of_data_matrix m =
  let (all_as :: [a]) = range fullRange in
  Matrix (Array.array fullRange
          [ ((i, j), v)
          | (i, r) <- zip all_as (Data.Matrix.toLists m)
          , (j, v) <- zip all_as r])

matrix_inverse_precise :: (Ix' a, Num v, Fractional v, Eq v) => Matrix a a v -> Matrix a a v
matrix_inverse_precise x = trace "inverting" $ of_data_matrix . (\(Right res) -> res) . Data.Matrix.inverse . to_data_matrix $ x

to_hmatrix :: forall a . (Ix' a) => Matrix a a Rat -> Numeric.LinearAlgebra.HMatrix.Matrix Double
to_hmatrix (Matrix m) =
  Numeric.LinearAlgebra.Data.fromLists
   $ [
     [ fromRational . toRational $ (m ! (i, j))
       | j <- range fullRange
     ]
     | i <- range fullRange
     ]

of_hmatrix :: forall a . Ix' a => Numeric.LinearAlgebra.HMatrix.Matrix Double -> Matrix a a Rat
of_hmatrix m =
  let (all_as :: [a]) = range fullRange in
  Matrix (Array.array fullRange
          [ ((i, j), fromRational . toRational $ v)
          | (i, r) <- zip all_as (Numeric.LinearAlgebra.Data.toLists m)
          , (j, v) <- zip all_as r])

matrix_inverse_hmatrix :: (Ix' a) => Matrix a a Rat -> Matrix a a Rat
matrix_inverse_hmatrix x = of_hmatrix . Numeric.LinearAlgebra.HMatrix.inv . to_hmatrix $ x

to_eigen_matrix :: forall a . (Ix' a) => Matrix a a Rat -> Data.Eigen.Matrix.Matrix Double CDouble
to_eigen_matrix (Matrix m) =
  Data.Eigen.Matrix.fromList
   [
     [ fromRational . toRational $ (m ! (toEnum (i - 1), toEnum (j - 1)))
       | j <- range fullRange
     ]
     | i <- range fullRange
     ]

of_eigen_matrix :: forall a . Ix' a => Data.Eigen.Matrix.Matrix Double CDouble -> Matrix a a Rat
of_eigen_matrix m =
  let (all_as :: [a]) = range fullRange in
  Matrix (Array.array fullRange
          [ ((i, j), fromRational . toRational $ v)
          | (i, r) <- zip all_as (Data.Eigen.Matrix.toList m)
          , (j, v) <- zip all_as r])

matrix_inverse_eigen :: (Ix' a) => Matrix a a Rat -> Matrix a a Rat
matrix_inverse_eigen x = of_eigen_matrix . Data.Eigen.Matrix.inverse . to_eigen_matrix $ x

matrix_inverse_verified inverse x =
  let y = inverse x in
    if (matrix_mult x y == Matrix (f_array (\(x, y) -> if x == y then 1 else 0))) then y else error $ "matrix inverse broken: " ++ "\n" ++ show x ++"\n" ++ show y

-- x = m * x + x_0
-- solve_equation :: (Ix' a, Num v, Fractional v, Eq v) => Matrix a a v -> Vector a v -> Vector a v
-- solve_equation a x0 = matrix_inverse (matrix_subtract a matrix_identity) `matrix_mult` matrix_negate x0

sparse_transpose :: (Ord a, Ord b) => Map a (Map b v) -> Map b (Map a v)
sparse_transpose m =
  fmap
  Map.fromList
  $ Map.fromListWith (++) (concatMap (\(a, m') -> map (\(b, v) -> (b, [(a, v)])) (Map.toList m')) (Map.toList m))

matrix_to_sparse :: (Ix' a, Ix' b, Ord b, Linear v, Eq v) => Matrix a b v -> [(a, (Map b v))]
matrix_to_sparse (Matrix m) =
  [ (a, Map.fromList
    [
      (b, x)
      | b <- range fullRange
      , let x = m ! (a,b)
      , x /= zero
    ])
    | a <- range fullRange
  ]

type RawProduct = Product

solvedRecipes
  :: GameConfig
     -> Map Product (Map RawProduct Rat, Map RecipeName Rat)
solvedRecipes gc =
  find_kernel_with_trace (/) (*) (recipesToMatrix gc)

currentSolvedRecipes = solvedRecipes current_game_config
currentRecipeMatrix = recipesToMatrix current_game_config

vector_lookup :: Ix' a => Vector a v -> a -> v
vector_lookup (Matrix x) a = x ! (a, ())

compute'_new :: GameConfig -> Product -> RawMaterialPressure
compute'_new gc =
  let recipes = fmap fst $ solvedRecipes gc in
   \product -> case Map.lookup product recipes of
     Nothing -> Map.empty
     Just m -> m

compute_recipe gc =
  let compute = compute'_new gc in
    \recipe ->
      mconcat'
        [ scale quantity (compute product)
        | (product, quantity) <- recipeProducts recipe
        ]

compute' gc = compute'_new gc

data SparseMatrix a b v = SparseMatrix [(a, Map b v)]

data MatrixError =
  Not_enough_equations
  | Contradicting_equations
  deriving Show

computeTotalCost :: Product -> RawMaterialPressure
computeTotalCost = compute' current_game_config

usability' Beacon = Usable
usability' GearWheel = Unusable
usability' IronPlate = Unusable
usability' IronOre = Unusable
usability' CopperOre = Unusable
usability' CopperCable = Unusable
usability' ElectronicCircuit = Unusable
usability' AdvancedCircuit = Unusable
usability' Plastic = Unusable
usability' Sulfur = Unusable
usability' SulfuricAcid = Unusable
usability' EngineUnit = Unusable
usability' CopperPlate = Unusable
usability' SteelPlate = Unusable
usability' ElectricMiningDrill = Usable
usability' Pipe = Usable
usability' Inserter = Usable
usability' ElectricFurnace = Usable
usability' SciencePack1 = Unusable
usability' SciencePack2 = Unusable
usability' SciencePack3 = Unusable
usability' AssemblingMachine1 = Usable
usability' AssemblingMachine2 = Usable
usability' AssemblingMachine3 = Usable
usability' TransportBelt = Usable
usability' StoneBrick = Unusable
usability' Lubricant = Unusable
usability' SciencePackProduction = Unusable
usability' SciencePackHighTech = Unusable
usability' CoalLiquefaction = Unusable
usability' ElectricEngineUnit = Unusable
usability' Stone = Unusable
usability' Coal = Unusable
usability' ProcessingUnit = Unusable
usability' SpeedModule = Usable
usability' SpeedModule2 = Usable
usability' SpeedModule3 = Usable
usability' EfficiencyModule = Usable
usability' EfficiencyModule2 = Usable
usability' EfficiencyModule3 = Usable
usability' ProductivityModule = Usable
usability' ProductivityModule2 = Usable
usability' ProductivityModule3 = Usable
usability' SolidFuel = Unusable
usability' ElectricalEnergy = Usable
usability' PetroleumGas = Unusable
usability' ResearchCoalLiquefaction = Unusable
usability' ResearchEndgame = Unusable
usability' ResearchRocketSilo = Unusable
usability' ResearchNuclearPower = Unusable
usability' ResearchLaserTurretDamage5 = Unusable
usability' LightOil = Unusable
usability' HeavyOil = Unusable
usability' CrudeOil = Unusable
usability' Steam = Unusable
usability' LaserTurret = Usable
usability' Battery = Unusable
usability' PiercingRoundMagazine = Usable
usability' FirearmRoundMagazine = Usable
usability' Grenade = Usable
usability' GunTurret = Usable
usability' SciencePackMilitary = Unusable
usability' SteelFurnace = Usable
usability' StoneFurnace = Usable
usability' LabBuilding = Usable
usability' SteamEngineBuilding = Usable
usability' BoilerBuilding = Usable
usability' OilRefinery = Usable
usability' ChemicalPlant = Usable
usability' SolarFacilityBuilding = Usable
usability' SolarPanel = Usable
usability' Accumulator = Usable
usability' Substation = Usable
usability' Roboport = Usable
usability' SteamTurbine = Usable
usability' NuclearReactor = Usable
usability' HeatExchanger = Usable
usability' HeatPipe = Usable
usability' Concrete = Usable
usability' RocketPart = Unusable
usability' RocketSiloBuilding = Unusable
usability' ControlModule = Unusable
usability' LightweightStructure = Unusable
usability' RocketFuel = Unusable
usability' SciencePackSpace = Unusable
usability' NuclearFacilityBuilding = Usable
usability' Satellite = Usable
usability' Radar = Usable
usability' ChemicalEnergy = Unusable
usability' BuriedIron = Unusable
usability' BuriedCopper = Unusable
usability' BuriedCoal = Unusable
usability' BuriedStone = Unusable
usability' Pollution = Unusable

usability recipeName =
  case recipeName of
    ProductRecipe product -> usability' product
    LiquefactionRecipe -> Unusable
    AdvancedOilProcessing -> Unusable
    UseAsFuelRecipe _product -> Unusable
    BoilerRecipe -> Usable

evaluateTotalCost  :: RawMaterialPressure -> Rat
evaluateTotalCost f = sum [ (estimate k * v) | (k, v) <- Map.toList f, v /= zero] where
--  estimate LightOil = 0.1
  estimate CrudeOil = 0.1
  estimate HeavyOil = 0.1
  estimate BuriedCoal = 1
  estimate BuriedIron = 1.5
  estimate BuriedCopper = 1
  estimate BuriedStone = 0.3
  estimate Pollution = 0.001
  estimate product = error $ "don't know how much this is worth: " ++ show product
--  estimate PetroleumGas = 0.1

subtract' a b = add a (minus b)

instance NFData Venue

data Assessment = Assessment {
  totalRawMaterials :: RawMaterialPressure,
  totalCapitalSeconds :: RawMaterialPressure,
  totalCapital :: RawMaterialPressure
} deriving (Generic, Show)

instance NFData Assessment

instance Linear Assessment where
  zero = Assessment {
    totalRawMaterials = zero,
    totalCapitalSeconds = zero,
    totalCapital = zero
    }
  a `add` b =
    Assessment
      {
        totalRawMaterials = add (totalRawMaterials a) (totalRawMaterials b),
        totalCapitalSeconds = add (totalCapitalSeconds a) (totalCapitalSeconds b),
        totalCapital = add (totalCapital a) (totalCapital b)
      }
  minus a =
    Assessment
      {
        totalRawMaterials = minus (totalRawMaterials a),
        totalCapitalSeconds = minus (totalCapitalSeconds a),
        totalCapital = minus (totalCapital a)
      }

instance VectorSpace Assessment where
  type Scalar Assessment = Rat
  scale x a =
    Assessment
      {
        totalRawMaterials = scale x (totalRawMaterials a),
        totalCapitalSeconds = scale x (totalCapitalSeconds a),
        totalCapital = scale x (totalCapital a)
      }

capitalUsePerExecution :: GameConfig -> Recipe -> RawMaterialPressure
capitalUsePerExecution gc recipe =
    let config = gc_configs gc recipe in
    let Time time = recipeTime recipe in
    let execution_time = time / (speedMultiplier config * baseSpeed gc (configVenue config)) in
    let facility_cost = capitalCost config in
    scale execution_time facility_cost

type CapitalUse = RawMaterialPressure

computeTotalCost_multi :: Map Product Rat -> RawMaterialPressure
computeTotalCost_multi = mconcat' . map (\(product, amount) -> scale amount (computeTotalCost product)) . Map.toList

assess_gc' :: GameConfig -> Map Product Rat -> (Map RecipeName CapitalUse, Assessment)
assess_gc' gc =
  let solved = solvedRecipes gc in
  let rbn = enabledRecipesByName (to_qualitative gc) in
  \demand -> runIdentity $ do
    (rawMaterials, executions) <- return $ dot_product_with scale demand solved
    capitalSeconds <- return $ Map.mapWithKey (\recipeName executions ->
        let Just recipe = Map.lookup recipeName rbn in
        scale executions (capitalUsePerExecution gc recipe)
      ) executions
    labResearchCapital <- return $ computeResearchCapital (gc_lab_researches_done gc) lab_speed_researches
    miningResearchCapital <- return $ computeResearchCapital (gc_mining_researches_done gc) mining_productivity_researches
    return $ (capitalSeconds, Assessment
      {
        totalRawMaterials = rawMaterials,
        totalCapitalSeconds = mconcat' (Map.elems capitalSeconds),
        totalCapital = computeTotalCost_multi (mconcat' [labResearchCapital, miningResearchCapital])
      })

assess_gc configs =
  let a = assess_gc' configs in
  \demand -> snd (a demand)

newtype Rat = Rat Double deriving (Eq, Ord, Generic, NFData, Linear, Num, Fractional, Real)

instance VectorSpace Rat where
  type Scalar Rat = Rat
  scale (Rat x) (Rat y) = Rat (x * y)

dot_product_with :: (Ord k, Linear a, Linear b, Linear c) => (a -> b -> c) -> Map k a -> Map k b -> c
dot_product_with f m1 m2 = foldr add zero $ Map.elems $
  Map.mergeWithKey
    (\_k a b -> Just (f a b))
    (\_ -> Map.empty)
    (\_ -> Map.empty)
    m1
    m2

instance Show Rat where
  show (Rat x) = Printf.printf "%.4f" (fromRational $ toRational x :: Double)

showModule SpeedModule = "s1"
showModule SpeedModule2 = "s2"
showModule SpeedModule3 = "s3"
showModule EfficiencyModule = "e1"
showModule EfficiencyModule2 = "e2"
showModule EfficiencyModule3 = "e3"
showModule ProductivityModule = "p1"
showModule ProductivityModule2 = "p2"
showModule ProductivityModule3 = "p3"
showModule x = error $ "not a module: " ++ show x

divv a b = if b == 0 then 1e10 else a / b

venueBuilding :: Venue -> [Product]
venueBuilding venue = case venue of
  Assembly2 -> [AssemblingMachine2]
  Assembly3 -> [AssemblingMachine3]
  Miner -> [ElectricMiningDrill]
  SmelterElectric -> [ElectricFurnace]
  SmelterBurner -> [SteelFurnace]
  Chemical -> [ChemicalPlant]
  Lab -> [LabBuilding]
  Boiler -> [BoilerBuilding]
  SteamEngine -> [SteamEngineBuilding]
  SolarFacility -> [SolarFacilityBuilding]
  NuclearFacility -> [NuclearFacilityBuilding]
  NoVenue -> []
  Refinery -> [OilRefinery]
  RocketSilo -> [RocketSiloBuilding]

capitalCost :: Config -> RawMaterialPressure
capitalCost config =
  computeTotalCost_multi (
    configModuleMaterials config
    `add`
    Map.fromListWith (+) (map (\p -> (p, 1)) (venueBuilding (configVenue config))))

partition_market :: (a -> (Rat, Rat)) -> [a] -> ([a], [a], [a], [a], [a])
partition_market evaluate l
  =
  ( p (\(gain, cost) -> gain >= 0 && cost <= 0) (\(_gain, cost) -> cost) -- free capital
  , p (\(gain, cost) -> gain >= 0 && cost <= 0) (\(gain, _cost) -> (-gain)) -- free resources
  , p (\(gain, cost) -> gain >= 0 && cost <= 0) (\(gain, cost) -> (gain * cost)) -- free both
  , p (\(gain, cost) -> gain > 0 && cost >= 0) (\(gain, cost) -> cost / gain) -- buy
  , p (\(gain, cost) -> gain < 0 && cost < 0) (\(gain, cost) -> gain / cost) -- sell
  ) where
  p predicate order = sortBy (comparing (order . evaluate)) $ filter (predicate . evaluate) l

allRecipeNames qgc = [recipeName recipe | (recipe, check) <- recipes, check qgc]

venueKind_of_venue venue = case venue of
  Assembly2 -> AssemblyVenueKind
  Assembly3 -> AssemblyVenueKind
  SmelterBurner -> SmelterVenueKind
  SmelterElectric -> SmelterVenueKind
  Chemical -> ChemicalVenueKind
  Miner -> MinerVenueKind
  Lab -> LabVenueKind
  Boiler -> BoilerVenueKind
  SteamEngine -> SteamEngineVenueKind
  Refinery -> RefineryVenueKind
  NoVenue -> NoVenueVenueKind
  SolarFacility -> GreenPowerVenueKind
  NuclearFacility -> GreenPowerVenueKind
  RocketSilo -> RocketSiloVenueKind

isVenueDefault venue =
  venue == currentDefaultVenue (venueKind_of_venue venue)

possibleSavings :: RawMaterialPressure -> GameConfig -> [(Change, GameConfig, Assessment)]
possibleSavings demand gc =
  let assessment_base = assess_gc gc demand in
  let
   assess_diff gc' =
    let assessment_tip = assess_gc gc' demand in
    add assessment_tip (minus assessment_base)
  in
   (`using` parListChunk 10 rdeepseq) $ map (\(change, config) -> (change, config, assess_diff config)) (gc_alternatives gc)

showChange :: Change -> (String, String)
showChange (ProductChange product (venue, modules, beacon)) =
    let
     showVenue SmelterElectric = "+"
     showVenue SmelterBurner = "-"
     showVenue Assembly2 = "-"
     showVenue Assembly3 = "+"
     showVenue venue | isVenueDefault venue = "" | otherwise = "??"

     showBeacon Nothing = ""
     showBeacon (Just b) = "!" ++ showModule b
    in
    (show product, showVenue venue ++ (concatMap showModule modules) ++ showBeacon beacon)
showChange (Other a b) = (a, b)

data TableEntry = TableEntry {
  te_name :: (String, String),
  te_efficiency :: Time,
  te_saving :: RawMaterialPressure,
  te_cost :: RawMaterialPressure,
  te_saving_and_cost :: (Rat, Rat)
  }

evaluate (Time totalTime)
  name (Assessment {
              totalRawMaterials,
              totalCapitalSeconds,
              totalCapital}) =
  let te_saving = minus totalRawMaterials in
  let te_cost = scale (recip totalTime) totalCapitalSeconds `add` totalCapital in
  let efficiency = Time (totalTime * (evaluateTotalCost te_cost / evaluateTotalCost te_saving)) in
  let
   te_saving_and_cost =
    (evaluateTotalCost te_saving, evaluateTotalCost te_cost + installationCost )
  in
  (TableEntry { te_name = name, te_efficiency = efficiency, te_saving, te_cost, te_saving_and_cost })

possibleSavings' demand time =
  [ evaluate time (showChange change) assessment
  | (change, _gc, assessment) <- possibleSavings demand current_game_config
  ]

installationCost = 40000

desiredMaterials =
  [ (ResearchEndgame,  1)
  , (PiercingRoundMagazine, 20000)
  , (ProductivityModule3, 200)
  ]

lookup0 m k = case Map.lookup k m of
  Nothing -> zero
  Just x -> x

currentRecipeVenue recipe =
  let config = (gc_configs current_game_config recipe) in
  configVenue config

currentEffectiveExecutionTime recipeName =
  let recipe = (enabledRecipesByName (to_qualitative current_game_config) !!! recipeName) in
  unTime (recipeTime recipe)
     / (speedMultiplier (gc_configs current_game_config recipe) * baseSpeed current_game_config (currentRecipeVenue recipe))

showHours (Time t) =
  show (t / 3600) ++ "h"

format_material_pressure = show . evaluateTotalCost

rCols =
  [ ("Efficiency", (showHours . te_efficiency))
  , ("Name", (fst . te_name))
  , ("Mod", (snd . te_name))
  , ("Gain", (format_material_pressure . te_saving))
  , ("Cost", (format_material_pressure . te_cost))
  ]

pad n l = replicate (n - length l) ' ' ++ l

printTable :: [[String]] -> [String]
printTable =
  map concat . transpose . map (\col -> let maxl = maximum (map length col) in map (pad (maxl + 1)) col) . transpose

printTableG :: [a] -> [(String, (a -> String))] -> IO ()
printTableG l cols =
  let title = map fst cols in
  let showA row = map (($row) . snd) cols in
  mapM_ putStrLn $ printTable (title : map showA l)

printRs l = printTableG l rCols

interestingProducts = []

show_percent n =
  Printf.printf "%.2f%%" (fromRational (toRational n) * (100 :: Double))

printDetailedCost cost =
    let overall = evaluateTotalCost cost in
    let
     r =
      map (\(k, v) ->
      let here = evaluateTotalCost (Map.singleton k v) in
      [show k, show v, show here, show_percent $ here / overall])
      (Map.toList cost) ++ [["Total", "_", show overall, "100.00%"]]
    in
    mapM_ putStrLn $ printTable r

print_config_details totalTime demand gc =
  let solved = solvedRecipes gc in
  let matrix = recipesToMatrix gc in
  let
   (total_cost_per_second, executions_per_second) =
    foldr add zero (
      map
        (\(product, amount) ->
            scale (recip $ unTime totalTime) (scale amount $ solved !!! product)
            ) desiredMaterials)
  in
  let
    negative_executions_per_second =
      filter
       ((<0) . snd)
        (Map.toList executions_per_second)
  in
  let
   effective_execution_time recipeName =
    let recipe = (enabledRecipesByName (to_qualitative gc)  !!! recipeName) in
    unTime (recipeTime recipe)
       / (speedMultiplier (gc_configs gc recipe) * baseSpeed gc (configVenue (gc_configs gc recipe)))
  in
  do
    mapM_ print negative_executions_per_second
    print "total factories:"
    let factories k = flip fmap (Map.lookup k executions_per_second) (* effective_execution_time k)
    let (capital_use_per_recipe, Assessment _raw total_capital_use total_capital_oneoff) = assess_gc' gc demand
    let scale_capital_use = scale (recip $ unTime totalTime)
    printTableG (sortBy (comparing factories) (allRecipeNames currentQGC)) $
      [ ("Name", show)
      , ("Factories", maybe "<none>" show . factories)
      , ("Price",
          \k ->
            case k of
              ProductRecipe product ->
                show $ evaluateTotalCost $ computeTotalCost product
              _ -> "<complex>")
      , ("Capital", (\k -> show $ evaluateTotalCost $ scale_capital_use (lookup0 capital_use_per_recipe k)))
      ] ++ flip map interestingProducts (\product -> (show product, (\k -> show $ lookup0 (matrix !!! k) product * (lookup0 executions_per_second k))))
    putStrLn "Total cost:"
    printDetailedCost (scale (unTime totalTime) total_cost_per_second)
    putStrLn "Total production capital:"
    print (evaluateTotalCost $ scale_capital_use total_capital_use)
    putStrLn "Total research capital:"
    print (evaluateTotalCost $ total_capital_oneoff)

report =
  let totalTime = Time (5 * 3600) in
  let
   demand =
    Map.fromList desiredMaterials
  in
  let savings = possibleSavings' demand totalTime in
   do
    let (free_capital, free_resources, free_both, buys, sells) = partition_market te_saving_and_cost savings
    print_config_details totalTime demand current_game_config
    putStrLn $ "Number of alternative configurations considered:" ++ show (length (gc_alternatives current_game_config))
    putStrLn "Free capital:"
    printRs (take 10 free_capital)
    putStrLn "Free resources:"
    printRs (take 10 free_resources)
    putStrLn "Free both:"
    printRs (take 10 free_both)
    putStrLn "Buys:"
    printRs (take 20 buys)
    putStrLn "Sells:"
    printRs (take 20 sells)

matrix_of_lists lists =
  Matrix (Array.array fullRange
          [ ((i, j), toRational v)
          | (i, r) <- zip (range fullRange) lists
          , (j, v) <- zip (range fullRange) r])

_identity_matrix :: (Ix' a) => Matrix a a Rat
_identity_matrix = Matrix (f_array (\(a,b) -> if a == b then 1 else 0))

currentDefaultVenue :: VenueKind -> Venue
currentDefaultVenue AssemblyVenueKind = Assembly2
currentDefaultVenue SmelterVenueKind = SmelterElectric
currentDefaultVenue GreenPowerVenueKind = NuclearFacility
currentDefaultVenue venueKind = case venuesByKind venueKind of
  [ venue ] -> venue
  _ -> error "ambiguous venue"

trivial recipe =
  (currentDefaultVenue (recipeVenueKind recipe), [], Nothing)

(!!!) :: (Ord k, HasCallStack) => Map k v -> k -> v
(!!!) m x = m Map.! x

parseSpec :: String -> (Maybe Char, [Product], Maybe Product)
parseSpec ('+' : rest) = (\(m,b) -> (Just '+', m, b)) (parseSpec' rest)
parseSpec ('-' : rest) = (\(m,b) -> (Just '-', m, b)) (parseSpec' rest)
parseSpec rest = (\(m,b) -> (Nothing, m, b)) (parseSpec' rest)
parseSpec' ('!' : rest) = ([], parseBeaconSpec rest)
parseSpec' (c1 : c2 : rest) = (parseModule c1 c2 : ms, b) where
  (ms, b) = parseSpec' rest
parseSpec' [] = ([], Nothing)
parseSpec' s = error $ "syntax error when parsing modules configuration " ++ show s
parseBeaconSpec [c1, c2] = Just (parseModule c1 c2)
parseBeaconSpec [] = Nothing
parseBeaconSpec s = error $ "syntax error when parsing beacon configuration " ++ show s
parseModule 'p' '3' = ProductivityModule3
parseModule 'p' '2' = ProductivityModule2
parseModule 'p' '1' = ProductivityModule
parseModule 's' '3' = SpeedModule3
parseModule 's' '2' = SpeedModule2
parseModule 's' '1' = SpeedModule
parseModule 'e' '3' = EfficiencyModule3
parseModule 'e' '2' = EfficiencyModule2
parseModule 'e' '1' = EfficiencyModule
parseModule c1 c2 = error $ "unrecognized module: " ++ [c1, c2]

venueByChar :: VenueKind -> Maybe Char -> Venue
venueByChar SmelterVenueKind (Just '+') = SmelterElectric
venueByChar SmelterVenueKind (Just '-') = SmelterBurner
venueByChar SmelterVenueKind Nothing = error "ambiguous smelter"
venueByChar AssemblyVenueKind (Just '-') = Assembly2
venueByChar AssemblyVenueKind (Just '+') = Assembly3
venueByChar AssemblyVenueKind Nothing = error "ambiguous assembly"
venueByChar GreenPowerVenueKind (Just '-') = SolarFacility
venueByChar GreenPowerVenueKind (Just '+') = NuclearFacility
venueByChar GreenPowerVenueKind Nothing = error "ambiguous green power"
venueByChar kind Nothing = currentDefaultVenue kind
venueByChar kind (Just c) = error $ "venueByChar weird invocation" ++ show (kind, c)

currentSpecs =
 let
  products =
   init[
    --(ProcessingUnit,"+s2p3p3p3"),
    --(GearWheel,"+p2p2p2s1"),
    --(Plastic,"p2p2e1"),
    (SulfuricAcid,"p2p2p2"),
    --(AdvancedCircuit,"+e1e1p1p1"),
    --(ResearchNuclearPower,"p2p2"),
    --(ResearchRocketSilo,"p2p2"),
    --(EngineUnit,"-e1e1"),
    (SciencePack1,"-e1e1"),
    (SciencePack2,"-e1e1"),
    (LightOil,"e1e1e1"),
    --(CopperCable,"+e1e1e1p1"),
    --(PetroleumGas,"e1e1e1"),
    (IronOre,"e1e1e1"),
    (CopperOre,"e1e1e1"),
    (Coal,"e1e1e1"),
    (PiercingRoundMagazine,"-e1e1"),
    --(SciencePackHighTech,"+p3p3p3p3"),
    --(ElectronicCircuit,"+p3p3p3s2"),
    --(SciencePack3,"+p3p3p3s2")

    (ProcessingUnit,"+p3p3p3p3!s3"),
    (SciencePack3,"+p3p3p3p3!s3"),
    (SciencePackHighTech,"+p3p3p3p3!s3"),
    --(CopperCable,"+p2p2p2p2!s2"),
    (ElectronicCircuit,"+p3p3p3p3!s2"),
    (ResearchEndgame,"p3p3!s2"),
    (GearWheel,"+p3p3p3p3!s3"),

    (AdvancedCircuit,"+p3p3p3p3!s3"),
    (SciencePackMilitary,"+p3p3p3p3!s2"),
    (Plastic,"p3p3p3!s3"),

    (CopperPlate,"+p2p2!s2"),
    (IronPlate,"+p2p2!s2"),
    (SteelPlate,"+p2p2!s2"),
    (CopperCable, "+p3p3p3p3!s3"),
    (PetroleumGas, "p2p2p2!s2"),

    (RocketPart, "p3p3p3p3"),
    (ControlModule, "+p3p3p3p3!s3"),
    (RocketFuel, "p3p3p3p3!s3"),
    (LightweightStructure, "+p3p3p3p3!s3"),

    -- make sure:
    (SciencePackProduction,"+p3p3p3p3!s2"),

    (EngineUnit,"+p3p3p3p3!s2"),
   undefined]
 in
 map (\(p,s) -> (ProductRecipe p, s)) products ++ [
   (AdvancedOilProcessing, "p3p3p3!s3")
 ]

currentModules :: RecipeName -> PreConfig
currentModules = runIdentity $ do
  specs <- return $ Map.fromListWith (error "spec for the same product given twice") (map (second parseSpec) currentSpecs)
  let byName = Map.fromList (map (\(r, _) -> (recipeName r, r)) recipes)
  return $ \recipeName ->
    let recipe = byName !!! recipeName in
    case Map.lookup recipeName specs of
      Nothing ->
        trivial recipe
      Just (venue_char, modules, beacon) ->
        let venue = venueByChar (recipeVenueKind recipe) venue_char in
        (venue, modules, beacon)

current_game_config :: GameConfig
current_game_config =
 GameConfig {
  gc_lab_researches_done = 4,
  gc_mining_researches_done = 6,
  gc_liquefaction = Liquefaction_disabled,
  gc_recipe_configs = currentModules
  }

currentQGC = to_qualitative current_game_config

--main = print $ computeTotalCost SciencePack3
main = report
--main =
--  print $ solvedRecipes currentConfig !!! SulfuricAcid
-- main = mapM_ print $ possibleSavings'''
{-
main = do
 flip mapM_ (range fullRange :: [Product]) $ \product -> do
  let configs =
        (\p -> if p == GearWheel
          then Config {configSpeedBonus = 0, configProductivityBonus = 0, configEnergyBonus = 0}
          else currentConfig p)
  let x = compute' configs GearWheel
  let y = compute' configs GearWheel
  print $ x
  print $ y
  print $ x == y
-}

max_stars_repo_name: Rotsor/factorio-module-selector

max_issues_repo_name: Rotsor/factorio-module-selector

max_forks_repo_name: Rotsor/factorio-module-selector


text:
module Pendulum.View (
  drawPendulum
  ) where

import qualified Data.Vector.Sized as V
import qualified Data.Vector.Storable as VS
import           GHC.TypeLits
import           Graphics.Gloss
import           Numeric.Hamilton
import           Numeric.LinearAlgebra.Static
import           Pendulum.Model (Pendulum(..))

-- | 振り子の長さの描画スケール
lengthScale :: Float
lengthScale = 100.0

-- | 振り子のおもりの描画スケール
massScale :: Float
massScale = 10.0

r2list :: KnownNat n => R n -> [Double]
r2list = VS.toList . extract

drawPendulum :: KnownNat n => Pendulum n -> Picture
drawPendulum pd =
  let ls  = map ((*lengthScale) . realToFrac) (V.toList $ _lengthes pd)
      ms  = map ((*massScale) . realToFrac) (V.toList $ _masses pd)
      ps  = map realToFrac (r2list (cfgPositions (_config pd)))
      ths = map (subtract (pi / 2.0)) ps
      xs = scanl (+) 0 (zipWith (*) ls (map cos ths))
      ys = scanl (+) 0 (zipWith (*) ls (map sin ths))
      lines   = line (zip xs ys)
      circles = zipWith3 (translate) (tail xs) (tail ys) (map circleSolid ms)
   in pictures (lines:circles)


text:
{-# OPTIONS_GHC -fno-warn-orphans #-}

module TinyLang.Prelude
    ( module Export
    , module TinyLang.Prelude
    ) where

-- base
------------------------------
import           Control.Applicative       as Export
import           Control.Exception         as Export
import           Control.Monad             as Export
import           Control.Monad.IO.Class    as Export
import           Data.Bifunctor            as Export
import           Data.Bool                 as Export
import           Data.Char                 as Export
import           Data.Coerce               as Export
import           Data.Complex              as Export
import           Data.Either               as Export
import           Data.Foldable             as Export
import           Data.Function             as Export
import           Data.Functor              as Export
import           Data.Functor.Compose      as Export
import           Data.Functor.Identity     as Export
import           Data.List                 as Export
import           Data.Maybe                as Export
import           Data.Monoid               as Export hiding (First (..),
                                                      Last (..))
import           Data.Ord                  as Export
import           Data.Proxy                as Export
import           Data.Ratio                as Export
import           Data.Semigroup            as Export
import           Data.String               as Export
import           Data.Traversable          as Export
import           Data.Tuple                as Export
import           Data.Void                 as Export
import           Debug.Trace               as Export
import           GHC.Exts                  as Export (groupWith, sortWith)
import           GHC.Generics              as Export (Generic, Generic1)
import           GHC.Stack                 as Export
import           Numeric                   as Export
import           System.IO                 as Export
import           Text.Read                 as Export (Read (..), readEither,
                                                      readMaybe)

-- transformers
------------------------------
import           Control.Monad.Trans.Maybe as Export

-- mtl
------------------------------
import           Control.Monad.Except      as Export
import           Control.Monad.Reader      as Export
import           Control.Monad.State       as Export

-- containers
------------------------------
import           Data.IntMap.Strict        as Export (IntMap)
import           Data.IntSet               as Export (IntSet)
import           Data.Map.Strict           as Export (Map)
import           Data.Set                  as Export (Set)

-- hashable
------------------------------
import           Data.Hashable             as Export

-- unordered-containers
------------------------------
import           Data.HashMap.Strict       as Export (HashMap)
import           Data.HashSet              as Export (HashSet)

import qualified Data.IntMap.Strict        as IntMap

import           Test.QuickCheck.Property

-- vector
------------------------------

import           Data.Vector               as Export (Vector)

infixr 9 .*
infixr 2 ?

newtype PairT b f a = PairT
    { unPairT :: f (b, a)
    }

instance Functor f => Functor (PairT b f) where
    fmap f (PairT p) = PairT $ fmap (fmap f) p

instance Hashable a => Hashable (IntMap a) where
    hashWithSalt salt = hashWithSalt salt . IntMap.toList

-- Feels so weird not to have it by default.
instance Testable (Either String ()) where
    property = property . \case
        Left err -> failed { reason = err }
        Right () -> succeeded

(.*) :: (c -> d) -> (a -> b -> c) -> a -> b -> d
(.*) = (.) . (.)
{-# INLINABLE (.*) #-}

(?) :: Alternative f => Bool -> a -> f a
(?) b x = x <$ guard b
{-# INLINABLE (?) #-}

visitExtract :: (Functor t, Functor f) => (t a -> a) -> (a -> f b) -> t a -> f (t b)
visitExtract ext f a = (<$ a) <$> f (ext a)
{-# INLINABLE visitExtract #-}

firstF :: Functor f => (a -> f c) -> (a, b) -> f (c, b)
firstF f (x, y) = flip (,) y <$> f x
{-# INLINABLE firstF #-}

foldlM' :: (Foldable t, Monad m) => (b -> a -> m b) -> b -> t a -> m b
foldlM' f z xs = foldr step return xs z where
  step x r a = f a x >>= (r $!)
{-# INLINABLE foldlM' #-}

foldMapA :: (Foldable t, Applicative f, Monoid b) => (a -> f b) -> t a -> f b
foldMapA f = foldr (\x r -> mappend <$> f x <*> r) (pure mempty)
{-# INLINABLE foldMapA #-}

-- | Fold a monadic function over a 'Foldable'. The monadic version of 'foldMap'.
foldMapM :: (Foldable f, Monad m, Monoid b) => (a -> m b) -> f a -> m b
foldMapM f xs = foldr step return xs mempty where
    step x r z = f x >>= \y -> r $! z `mappend` y
{-# INLINABLE foldMapM #-}

asumMap :: (Foldable t, Alternative f) => (a -> f b) -> t a -> f b
asumMap f = foldr ((<|>) . f) empty
{-# INLINABLE asumMap #-}

mapMaybeA :: Applicative f => (a -> f (Maybe b)) -> [a] -> f [b]
mapMaybeA f = foldr (\x r -> maybe id (:) <$> f x <*> r) (pure empty)
{-# INLINABLE mapMaybeA #-}

fromOption :: (Foldable f, Applicative g) => g a -> f a -> g a
fromOption = foldr (const . pure)
{-# INLINABLE fromOption #-}

-- | This function generalizes 'eitherToMaybe', 'eitherToList',
-- 'listToMaybe' and other such functions.
reoption :: (Foldable f, Alternative g) => f a -> g a
reoption = fromOption empty
{-# INLINABLE reoption #-}

forBind :: (Monad m, Traversable m, Applicative f) => m a -> (a -> f (m b)) -> f (m b)
forBind a f = join <$> traverse f a
{-# INLINABLE forBind #-}

max_stars_repo_name: effectfully/tiny-lang

max_issues_repo_name: effectfully/tiny-lang

max_forks_repo_name: jakzale/tiny-lang


text:
{-# OPTIONS_GHC  -fno-warn-unused-binds -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-missing-signatures #-}
{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, UndecidableInstances, FlexibleContexts, TypeSynonymInstances #-}

---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------
-- |
-- | Module : Approximate median
-- | Creator: Xiao Ling
-- | Created: 12/08/2015
-- | TODO   : test standard deviation of alpha, beta, and final version
-- |
---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------

module ApproxMedian (

      MaxSize
    , Store

    , aMedian
    , aMedian'

  ) where

import Control.Monad.Trans
import Control.Monad.State
import Control.Monad.Reader
import Control.Monad.Identity

import Data.Random
import Data.Conduit
import Data.Sequence
import Data.Foldable (toList)
import qualified Data.Conduit.List as Cl

import Core
import Statistics

{-----------------------------------------------------------------------------
  Types
------------------------------------------------------------------------------}

type MaxSize = Int
type Store a = (Seq a, MaxSize)

-- * Construct an empty store given max size `s`
store :: MaxSize -> Store a
store = (,) empty

{-----------------------------------------------------------------------------
  Median II
------------------------------------------------------------------------------}

-- * determine the `e`ps -approx median of a list `xs` with confidence `d`elta.
aMedian :: (Ord a, Floating a) => Eps -> Delta -> [a] -> IO a
aMedian e d xs = eval $ Cl.sourceList xs $$ toMedian $ store s
  where eval p = fmap (median . toList) . fmap fst . flip runRVar StdRandom $ evalStateT p 0
        s      = round $ 7/(e^2) * log (2/d)

-- * `tick` up a counter for each item `a` seen
-- * and with probability min(1, s/i) put item `a` into store `t`
toMedian :: Store a -> Sink a (StateT Counter RVar) (Store a)
toMedian t = Cl.foldM (\t a -> tick >> t $|> a) t

-- * count items seen so far
tick :: Enum s => MonadState s m => m ()
tick = modify succ

-- * with probability `s/i` uniformly select an item from the store `t` and replace it with `a`
-- * mnemonic: `|>` is insertion for Sequences, and `$f` means do f probabilistically
($|>) :: Store a -> a -> StateT Counter RVar (Store a)
($|>) t@(as,s) a = do

  i <- get

  let p = s ./ i
  if p >= 1 then return (as |> a, s) else do

      h <- lift . toss $ coin p

      if isHead h then do
          x <- lift (uniform 0 s :: RVar Int)
          let as' = update x a as
          return (as',s)
      else return t

-- * Sledghammer division
(./) :: (Fractional a, Read a, Show s, Show i) => s -> i -> a
(./) s i = (read . show $ s)/(read . show $ i)

{-----------------------------------------------------------------------------
  Median II
------------------------------------------------------------------------------}

-- * determine the `e`ps - approx median of a list `xs` with confidence `d`elta.
aMedian' :: (Ord a, Floating a) => Eps -> Delta -> [a] -> IO a
aMedian' e d xs = eval $ Cl.sourceList xs $$ toMedian'
  where eval p = fmap (median . toList) . fmap fst $ runRVar (evalStateT (runStateT p 0) (store s)) StdRandom
        s      = round $ 7/(e^2) * log (2/d)

-- * `tick` up a counter `i` of items seen so far
-- * and with uniform probability put the item `a` into a store `as` of max size `s`
toMedian' :: Sink a (StateT Counter (StateT (Store a) RVar)) (Seq a)
toMedian' = tick' >> do
  ma <- await
  case ma of
    Nothing -> (lift . lift $ get) >>= return . fst
    Just a  -> (lift . place $ a ) >>  toMedian'

tick' :: (Enum s, MonadTrans t, MonadState s m) => t m ()
tick' = lift . modify $ succ

-- * with uniform probability `s/i` select an item from the
-- * store and re`place` it with `a`
place :: a -> StateT Counter (StateT (Store a) RVar) ()
place a = do

  i        <- get
  t@(as,s) <- lift get

  let p = s ./ i

  if p >= 1 then lift $ put (as |> a, s) else do

      h <- lift . lift . toss $ coin p
      case isHead h of
        False -> lift $ put (as,s)
        _     -> do
          x <- lift $ lift (uniform 0 s :: RVar Int)
          let as' = update x a as
          lift $ put (as',s)

{-----------------------------------------------------------------------------
  Awaiting SO response for this stuff
------------------------------------------------------------------------------}

instance MonadState s m => MonadState s (RVarT m) where
  get   = lift get
  put   = lift . put
  state = lift . state

foo :: (MonadState s m, MonadRandom m) => s -> RVarT m ()
foo s = do
  x <- (uniformT 0 1 :: RVarT m Double)
  if x < 0.5 then put s else return ()

foo' :: (MonadState String m, MonadRandom m) => RVarT m ()
foo' = foo "hello"

--foo'' :: (MonadState String m, MonadRandom m) => m ()
--foo''= runRVarT foo' StdRandom

bar :: (MonadTrans t1, MonadState t m, MonadState s (t1 m)) => s -> RVarT (t1 m) ()
bar s = do
  x <- (uniformT 0 1 :: RVarT m Double)
  a <- lift . lift $ get
  if x < 0.5 then put s else return ()

-- * desired behavior with state and ranndomness

baz :: (MonadState s m, MonadReader t m) => m ()
baz = get >> ask >> return ()

baz1 :: MonadReader t m => m ((),Int)
baz1 = runStateT baz 0

baz1' :: Identity ((),Int)
baz1' = runReaderT baz1 (1,"hello")

baz2 :: MonadState s m => m ()
baz2 = runReaderT baz "hello"

baz2' :: Identity ((),[Int])
baz2' = runStateT baz2 [1..10]


text:
{-# LANGUAGE DeriveGeneric, DeriveAnyClass #-}
{-# LANGUAGE FlexibleContexts #-}
module Neural.Layer
  ( Layer(..)
  , newLayer
  , randomLayer
  , feed
  )
where

import           System.Random
import           Data.Random.Normal
import           Control.Monad
import           Control.Monad.State
import           Numeric.LinearAlgebra
import           GHC.Generics                   ( Generic )
import           Control.DeepSeq

import           Neural.Activation              ( ActivationFunction(..) )

-- | A layer in a neural network.
data Layer a = Layer
    { weights :: Matrix a    -- Weight matrix
    , biases  :: Vector a    -- Bias vector
    } deriving (Show, Generic, NFData)

-- | Construct a new layer from a weigth matrix and a bias vector.
newLayer :: Matrix a -> Vector a -> Layer a
newLayer ws bs = Layer { weights = ws, biases = bs }

-- | Construct a new random layer given its,
-- size and the size of the previous layer.
-- Each bias and weight is given a random value
-- normally distributed with mean 0 and standard deviation 1.
randomLayer
  :: (RandomGen g, Random a, Element a, Floating a)
  => Int         -- ^ The size of the previous layer
  -> Int         -- ^ The size of the layer
  -> State g (Layer a)
randomLayer m n = do
  _biases  <- replicateM n $ state normal
  _weights <- replicateM n $ replicateM m $ state normal
  return $ Layer (fromLists _weights) (fromList _biases)

-- | Calculate the activation and weighted input of a layer
feed
  :: (Numeric a, Num (Vector a))
  => ActivationFunction a  -- ^ The activation function
  -> Vector a              -- ^ The layer input
  -> Layer a               -- ^ The layer
  -> (Vector a, Vector a)  -- ^ The layer output
feed _σ x (Layer w b) = (σ y, y)
 where
  y = w #> x + b
  σ = cmap _σ

max_stars_repo_name: cornelius-sevald/nnhd

max_issues_repo_name: cornelius-sevald/nnhd

max_forks_repo_name: cornelius-sevald/nnhd