text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE ExplicitNamespaces  #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE FlexibleInstances   #-}
{-# LANGUAGE IncoherentInstances #-}
{-# LANGUAGE InstanceSigs        #-}
{-# LANGUAGE OverloadedStrings   #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE Strict              #-}
{-# LANGUAGE TypeFamilies        #-}
{-# LANGUAGE TypeOperators       #-}
{-# LANGUAGE Unsafe              #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}

module Releaser.Build
    ( buildBORLTable
    , buildBORLGrenade
    , buildSim
    , borlSettings
    , scaleAlg
    , nnConfig
    , netInp
    , modelBuilder
    , actionConfig
    , experimentName
    , load
    , save
    , mInverse
    , databaseSetting
    , expSetting
    , mkMiniPrettyPrintElems
    ) where

import           Control.Arrow                     (second)
import           Control.DeepSeq
import           Control.Lens
import           Control.Monad
import           Control.Monad.IO.Class
import qualified Data.ByteString                   as B
import           Data.Constraint                   (Dict (..))
import           Data.Default
import           Data.Int                          (Int64)
import           Data.List                         (find, genericLength)
import qualified Data.Map                          as M
import           Data.Maybe                        (fromMaybe, isJust)
import qualified Data.Proxy                        as Px
import           Data.Reflection                   (reifyNat)
import           Data.Serialize                    as S
import           Data.Serialize.Text               ()
import qualified Data.Text                         as T
import qualified Data.Text                         as T
import qualified Data.Text.Encoding                as E
import qualified Data.Vector.Storable              as V
import           GHC.Exts                          (IsList (..))
import           GHC.TypeLits                      (KnownNat)
import           Grenade
import           Network.HostName
import           Prelude
import           Statistics.Distribution
import           Statistics.Distribution.Uniform
import           System.Directory
import           System.Environment                (getArgs)
import           System.IO.Unsafe                  (unsafePerformIO)
import           Unsafe.Coerce                     (unsafeCoerce)

-- ANN modules
import           Grenade

import           Experimenter                      hiding (sum)
import           ML.BORL                           as B hiding (actionFilter,
                                                         featureExtractor)
import qualified ML.BORL                           as B
import           SimSim                            hiding (productTypes)

import           Releaser.Costs.Type
import           Releaser.Decay.Type
import           Releaser.FeatureExtractor.Type
import           Releaser.Release.ReleasePlt
import           Releaser.Routing.Type
import           Releaser.SettingsAction
import           Releaser.SettingsActionFilter
import           Releaser.SettingsConfigParameters
import           Releaser.SettingsCosts
import           Releaser.SettingsDecay
import           Releaser.SettingsDemand
import           Releaser.SettingsFeatureExtractor
import           Releaser.SettingsPeriod
import           Releaser.SettingsReward
import           Releaser.SettingsRouting
import           Releaser.Type
import           Releaser.Util

import           Debug.Trace

buildSim :: IO SimSim
buildSim =
  newSimSimIO
    (configRoutingRoutes routing)
    -- procTimesConst
    (configProcTimes procTimes)
    periodLength
    -- (releaseBIL $ M.fromList [(Product 1, 1), (Product 2, 1)])
    -- (releaseBIL $ M.fromList [(Product 1, 2), (Product 2, 2)])
    -- (releaseBIL $ M.fromList [(Product 1, 3), (Product 2, 3)])
    -- (releaseBIL $ M.fromLisqt [(Product 1, 4), (Product 2, 4), (Product 3, 4), (Product 4, 4), (Product 5, 4), (Product 6, 4)])
    -- (releaseBIL $ M.fromList (map (\pt -> (pt, 2)) productTypes))
    -- releaseImmediate
    (mkReleasePLT initialPLTS)
    dispatchFirstComeFirstServe
    shipOnDueDate

initialPLTS :: M.Map ProductType Time
initialPLTS = M.fromList $ zip productTypes (replicate len 1)
  where
    len
      | bnNbn = 2
      | otherwise = length productTypes

-- testDemand :: IO ()
-- testDemand = do
--   let nr = 1000
--   g <- createSystemRandom
--   sim <- buildSim
--   xs <- replicateM nr (generateOrders sim)
--   let len = fromIntegral $ length (concat xs)
--   putStr "Avg order slack time: "
--   print $ timeToDouble (sum $ map orderSlackTime (concat xs)) / len
--   putStr "Avg order arrival date: "
--   print $ timeToDouble (sum $ map arrivalDate (concat xs)) / len
--   putStr "Avg number of order per period: "
--   print $ fromIntegral (sum (map length xs)) / fromIntegral nr
--   putStr "Avg order due date: "
--   print $ timeToDouble (sum $ map dueDate (concat xs)) / len
--   putStr "Avg number of order per product type"
--   print $ map length $ groupBy ((==) `on` productType) $ sortBy (compare `on` productType) (concat xs)
--   print $ map (productType . head) $ groupBy ((==) `on` productType) $ sortBy (compare `on` productType) (concat xs)

------------------------------------------------------------
--------------------------- BORL ---------------------------
------------------------------------------------------------

instance Serialize Release where
  put (Release _ n) = S.put $ T.unpack n
  get = do
    n <- T.pack <$> S.get
    let fun | n == pltReleaseName = mkReleasePLT initialPLTS
            | n ==  uniqueReleaseName releaseImmediate = releaseImmediate
            | T.isPrefixOf bilName n = releaseBIL $ M.fromList bilArgs
              where ~bilArgs = read (T.unpack $ T.drop (T.length bilName) n)
                    ~bilName = T.takeWhile (/= '[') $ uniqueReleaseName (releaseBIL mempty)
    return fun

instance Show St where
  show st = show (extractFeatures False st)

netInp :: St -> V.Vector Double
netInp = extractionToList . extractFeatures True

mInverse :: BORL St Act -> NetInputWoAction -> Maybe (Either String St)
mInverse borl = return . Left . show . fromListToExtraction (borl ^. s) (featureExtractor True)

netInpTbl :: St -> V.Vector Double
netInpTbl st = case extractFeatures False st of
  Extraction plts op que _ fgi shipped _ -> V.fromList $ plts ++ map reduce (concat $ op ++ map (map (fromIntegral . ceiling . (/9))) (concat que) ++ fgi ++ shipped)
  where
    reduce x = 7 * fromIntegral (ceiling (x / 7))

netInpTblBinary :: St -> [Double]
netInpTblBinary st = case extractFeatures False st of
  Extraction plts op que _ fgi shipped _ -> plts ++ map reduce (concat op) ++ map (fromIntegral . ceiling . (/9)) (concat (concat que)) ++ map reduce (concat $ fgi ++ shipped)
  where
    reduce x | x == 0 = x
             | otherwise = 1

modelBuilder :: NrFeatures -> (NrRows, NrCols) -> IO SpecConcreteNetwork
modelBuilder nrInp (nrRows, nrCols) =
  buildModelWith (NetworkInitSettings UniformInit BLAS Nothing) (DynamicBuildSetup False) $
  inputLayer1D nrInp >>
  fullyConnected (round $ 1.75*fromIntegral nrInp) >> leakyRelu >>
  fullyConnected (round $ 1.5* fromIntegral nrInp) >> leakyRelu >>
  fullyConnected lenOut >> reshape (nrRows, nrCols, 1) -- >> tanhLayer -- leakyTanhLayer 0.98
  where
    lenOut = nrRows * nrCols

-- SpecFullyConnected 180 540 :=> SpecRelu (540,1,1) :=> SpecFullyConnected 540 180 :=> SpecRelu (180,1,1) :=> SpecFullyConnected 180 24 :=> SpecRelu (24,1,1) :=> SpecFullyConnected 24 12 :=> SpecReshape (12,1,1) (6,2,1) :=> SpecLeakyTanh 0.98 (6,2,1) :=> SpecNNil2D 6x2

--  *GOOD RESULTS* with DDS=12 and
-- SpecFullyConnected 45 90 :=> SpecRelu (90,1,1) :=> SpecFullyConnected 90 45 :=> SpecRelu (45,1,1) :=> SpecFullyConnected 45 6 :=> SpecRelu (6,1,1) :=> SpecFullyConnected 6 3 :=> SpecTanh (3,1,1) :=> SpecNNil1D 3

mkInitSt :: SimSim -> [Order] -> (AgentType -> IO St, St -> [V.Vector Bool])
mkInitSt sim startOrds =
  let initSt = St sim startOrds rewardFunction initialPLTS
      actFilter = mkConfig actionFilter actionFilterConfig
  in (return . const initSt, actFilter)

buildBORLTable :: IO (BORL St Act)
buildBORLTable = do
  sim <- buildSim
  startOrds <- liftIO $ generateOrders sim
  let (initSt, actFilter) = mkInitSt sim startOrds
  mkUnichainTabular alg initSt netInpTbl -- netInpTblBinary
    action actFilter borlParams (configDecay decay) borlSettings (Just initVals)

buildBORLGrenade :: IO (BORL St Act)
buildBORLGrenade = do
  sim <- buildSim
  startOrds <- liftIO $ generateOrders sim
  let (initSt, actFilter) = mkInitSt sim startOrds
  st <- liftIO $ initSt MainAgent
  flipObjective . setPrettyPrintElems <$>
    mkUnichainGrenadeCombinedNet alg initSt netInp action actFilter borlParams (configDecay decay) modelBuilder nnConfig borlSettings (Just initVals)
  -- flipObjective . setPrettyPrintElems <$> mkUnichainGrenade alg initSt netInp action actFilter borlParams (configDecay decay) (modelBuilder st) nnConfig borlSettings (Just initVals)

setPrettyPrintElems :: BORL St Act -> BORL St Act
setPrettyPrintElems borl = setAllProxies (proxyNNConfig . prettyPrintElems) (ppElems borl) borl
  where ppElems borl = mkMiniPrettyPrintElems (borl ^. s)

copyFiles :: String -> ExperimentNumber -> RepetitionNumber -> Maybe ReplicationNumber -> IO ()
copyFiles pre expNr repetNr mRepliNr = do
  let dir = "results/" <> T.unpack (T.replace " " "_" experimentName) <> "/data/"
  createDirectoryIfMissing True dir
  mapM_
    (\fn -> copyIfFileExists fn (dir <> pre <> fn <> "_exp_" <> show expNr <> "_rep_" <> show repetNr <> maybe "" (\x -> "_repl_" <> show x) mRepliNr))
    ["reward", "stateValues", "episodeLength", "plts", "costs", "stateVAllStates", "stateWAllStates", "statePsiVAllStates", "statePsiWAllStates"]

copyIfFileExists :: FilePath -> FilePath -> IO ()
copyIfFileExists fn target = do
  exists <- doesFileExist fn
  when exists $ copyFileWithMetadata fn target

databaseSetting :: IO DatabaseSetting
databaseSetting = do
  hostName <- getHostName
  args <- getArgs
  let mHostArg = case find ((== "--host=") . take 7) args of
                   x@Just{} -> drop 7 <$> x
                   Nothing  -> drop 3 <$> find ((== "-h=") . take 3) args
  let getPsqlHost h
        | isJust mHostArg = maybe "" (E.encodeUtf8 .T.pack) mHostArg
        | h `elem` ["schnecki-zenbook", "schnecki-laptop"] = "192.168.1.110"
        | otherwise = "c437-pc147"
  putStrLn $ "Using DB-Host: " <> show (getPsqlHost hostName)
  return $ DatabaseSetting ("host=" <> getPsqlHost hostName <> " dbname=experimenter user=experimenter password=experimenter port=5432") 10

mkMiniPrettyPrintElems :: St -> [V.Vector Double]
mkMiniPrettyPrintElems st
  | length (head xs) /= length base' = error $ "wrong length in mkMiniPrettyPrintElems: " ++
                                show (length $ head xs) ++ " instead of " ++ show (length base') ++ ". E.g.: " ++ show (map (unscaleDouble scaleAlg (Just (scaleOrderMin, scaleOrderMax))) base') ++
                                "\nCurrent state: " ++ show (extractFeatures True st)

  | otherwise = map V.fromList $ concatMap (zipWith (++) plts . replicate (length plts) . map (scaleDouble scaleAlg (Just (scaleOrderMin, scaleOrderMax)))) xs
  where
    len = V.length $ extractionToList $ extractFeatures True st
    base' = drop lenPt (V.toList $ netInp st)
    lenPt | bnNbn = 2
          | otherwise = length productTypes

    plts :: [[Double]]
    plts = map (map (scaleDouble scaleAlg (Just (scalePltsMin, scalePltsMax))) . take lenPt) [[1, 3, 1, 1, 3, 1], [3, 5, 3, 3, 5, 3]]
    xs :: [[Double]]
    xs | len - length (head plts) == 22 = [xsSimple, xsSimple2]
       | len - length (head plts) == 21 = map init [xsSimple, xsSimple2]
       | len - length (head plts) == 18 = [xs19, xs19']
       | len - length (head plts) == 35 = [xsFull1, xsFull2]
       | len - length (head plts) == 30 = [xs30', xs30]
       | len - length (head plts) == 29 = [xs29', xs29]
       | len - length (head plts) == 50 = [init xs51]
       | len - length (head plts) == 51 = [xs51]
       | len - length (head plts) == 45 = [xs45]
       | len - length (head plts) == 44 = [xs44, xs44']
       | len - length (head plts) == 106 = [xs106, xs106']
       | len - length (head plts) == 156 = [xs156, xs156']
       | len - length (head plts) == 226 = [xs226]
       | len - length (head plts) == 23 = [xs23]
       | len - length (head plts) == 534 = [xs534]
       | len - length (head plts) == 138 = [xs138]
       | len - length (head plts) == 40 = [xs40]
       | trace ("len - length (head plts): " ++ show (len - length (head plts)))len - length (head plts) == 678 = [xs684]
       --  | otherwise = error ("No mkMiniPrettyPrintElems in Build.hs setup for length: " ++ show len ++ "\nCurrent state: " ++ show (extractFeatures True st))
       | otherwise = [replicate (len - length (head plts)) 0]
    xsSimple =              concat [concat [[ 0, 6, 8, 4, 4, 9, 9]], concat [ concat [[ 2]        ]], concat [[ 1]], concat [[ 0, 0, 0, 0, 0, 0]], concat [[ 0, 0, 0, 0, 0, 5, 4]]]
                            -- concat [concat [[ 0, 6, 8, 4, 4, 9, 9]], concat [ concat [[ 2],  [2],[2]  ]], concat [[ 1]], concat [[ 0, 0, 0, 0, 0, 0]], concat [[ 0, 0, 5, 4]]]
    xsSimple2 =             concat [concat [[ 0, 0, 7,10,10, 9,14]], concat [ concat [[23]         ]], concat [[ 1]], concat [[ 0, 0, 0, 0, 0, 0]], concat [[ 0, 0, 0, 3, 4, 9, 0]]
                                   -- concat [[ 0, 0, 7,10,10, 9,14]], concat [ concat [[23],  [2],[12]  ]], concat [[ 1]], concat [[ 0, 0, 0, 0, 0, 0]], concat [[ 0, 0, 9, 0]]
                                   ]
    xsFull1 = concat [ concat [[ 0, 6, 8, 4, 4, 9, 9]], concat  [concat [[ 2, 0, 0, 0, 0, 0, 0, 0]]], concat [[ 1, 0, 0, 0, 0, 0, 0, 0]], concat [[ 2, 0, 0, 0, 0, 0]], concat [[ 0, 0, 0, 0, 5, 4]]]
    xsFull2 = concat [ concat [[ 0, 0, 7, 10, 10, 9, 14]], concat  [concat [[ 2, 12, 9, 0, 0, 0, 0, 0]]], concat [[ 1, 0, 0, 0, 0, 0, 0, 0]], concat [[ 0, 0, 0, 0, 0, 0]], concat [[ 0, 0, 3, 4, 9, 0]]]
    xs51 = concat [ concat [[ 0, 0, 0, 0, 0, 0, 0,13]], concat [ concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], concat [[ 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0]]]
    xs45 = concat [ concat [[ 0, 0, 0, 0, 0, 0, 0,15]], concat [ concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]], concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0]]]
    xs30 = concat [concat [[ 0, 0, 7,10,10, 9,14]], concat [ concat [[ 0, 0, 0, 0, 0, 4, 6, 0, 0, 0, 0, 0]]], concat [[]], concat [[ 1, 0, 0, 0, 0, 0, 0]], concat [[  0, 0, 0, 0]]]
    xs30' = concat [concat [[ 0, 0, 7,10,10, 9,14]], concat [ concat [[ 0, 0, 0, 0, 0, 4, 6, 0, 0, 0, 0, 0]]], concat [[]], concat [[ 3, 2, 0, 0, 0, 0, 0]], concat [[ 3, 4, 9, 2]]]
    xs29 = concat [concat [[ 0, 0, 7,10,10, 9,14]], concat [ concat [[ 0, 0, 0, 0, 0, 4, 6, 0, 0, 0, 0, 0]]], concat [[]], concat [[ 1, 0, 0, 0, 0, 0]], concat [[  0, 0, 0, 0]]]
    xs29' = concat [concat [[ 0, 0, 7,10,10, 9,14]], concat [ concat [[ 0, 0, 0, 0, 0, 4, 6, 0, 0, 0, 0, 0]]], concat [[]], concat [[ 3, 2, 0, 0, 0, 0]], concat [[ 3, 4, 9, 2]]]
    xs19 = concat [concat [[ 0, 0, 7,10,10, 9,14]], concat [ concat [[ 4]]], concat [[]], concat [[ 1, 0, 0, 0, 0, 0]], concat [[  0, 0, 0, 0]]]
    xs19' = concat [concat [[ 0, 0, 7,10,10, 9,14]], concat [ concat [[ 6 ]]], concat [[]], concat [[ 3, 2, 0, 0, 0, 0]], concat [[ 3, 4, 9, 2]]]

    xs44 = concat [ concat [[ 0, 0, 12, 14, 15, 11, 4, 10, 12, 8, 3,13]], concat [ concat [[ 0, 0, 0, 0, 0, 4, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]],[], concat [[ 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], concat [[ 0, 0, 0, 2]]]
    xs44' = concat [ concat [[ 0, 0, 1, 4, 5, 1, 4, 10, 12, 15, 12,13]], concat [ concat [[ 0, 0, 0, 0, 0, 4, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]],[], concat [[ 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], concat [[ 0, 0, 0, 2]]]
    xs156 = concat [ concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8]],concat [ concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]],[], concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], concat [[ 0, 0, 0, 0],[ 0, 0, 0, 0]]]
    xs156' = concat [concat [[ 0, 0, 4, 8, 5, 6, 5, 8, 8, 2, 3, 1],[ 0, 0, 0, 0, 0, 0, 0, 4, 7, 3, 5, 2]],concat [ concat [[ 0, 0, 7, 3, 3, 7, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 3, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 7, 3, 3, 6, 5, 3, 2, 0, 0, 0, 0, 0, 0]]],[],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 1, 0, 0],[ 0, 0, 0, 5]]]

    xs106 = concat [ concat[[ 0, 0, 0, 0, 0, 0, 3],[ 0, 0, 0, 0, 0, 0, 7]],concat [ concat[[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]],[],concat [[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0],[ 0, 0, 0, 0]]]
    xs106' = concat[ concat [[ 0, 0, 3, 9, 4, 6, 1],[ 0, 0, 0, 0, 0, 0, 2]],concat [ concat [[ 0, 0, 0, 0, 0, 4, 5, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0,11, 2, 0]],concat [[ 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0]]],[],concat [[ 2, 0, 0, 0, 0, 0],[ 5, 1, 8, 0, 0, 0]],concat [[ 0, 0, 0, 0],[ 0, 0, 0, 0]]]

    xs684 :: [Double]
    xs684 = concat [ concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 2],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 3],[ 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, 1],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], concat[ concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]],[], concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0]]]
    xs226 :: [Double]
    xs226 = concat [ concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 5],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 3]], concat [ concat[[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]],[],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0],[ 0, 0, 0, 0]]]

    xs23 = concat [ concat [[ 0, 0, 0, 0, 0, 0, 6]],concat[concat[[ 0]],concat [[ 0]],concat [[ 0]],concat [[ 0]],concat [[ 0]],concat [[ 0]]],[],concat [[ 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0]]]
    xs534 = concat [ concat [[ 0, 0, 0, 0, 0, 0, 1],[ 0, 0, 0, 0, 0, 0, 2],[ 0, 0, 0, 0, 0, 0, 3],[ 0, 0, 0, 0, 0, 0, 1],[ 0, 0, 0, 0, 0, 0, 1],[ 0, 0, 0, 0, 0, 0, 0]], concat [concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]],[],concat [[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0]]]

    xs138 = concat [concat [[ 0, 0, 0, 2, 3, 1, 1],[ 0, 0, 0, 1, 2, 4, 2],[ 0, 0, 0, 1, 3, 2, 1],[ 0, 0, 0, 0, 1, 3, 3],[ 0, 0, 0, 2, 4, 2, 2],[ 0, 0, 0, 4, 2, 1, 4]],concat [concat [[ 0],[ 0],[ 0],[ 0],[ 0],[ 0]],concat [[ 0],[ 0],[ 0],[ 0],[ 0],[ 0]],concat [[ 0],[ 0],[ 0],[ 0],[ 0],[ 0]],concat [[ 0],[ 0],[ 0],[ 0],[ 0],[ 0]],concat [[ 0],[ 0],[ 0],[ 0],[ 0],[ 0]],concat [[ 0],[ 0],[ 0],[ 0],[ 0],[ 0]]],[],concat [[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0],[ 0, 0, 0, 0]]]

    xs40 = concat [concat [[ 0, 0, 0, 0, 0, 0, 4],[ 0, 0, 0, 0, 0, 0, 7]],concat [ concat [[ 0],[ 0]],concat [[ 0],[ 0]],concat [[ 0],[ 0]]],[],concat [[ 0, 0, 0, 0, 0, 0],[ 0, 0, 0, 0, 0, 0]],concat [[ 0, 0, 0, 0],[ 0, 0, 0, 0]]]

------------------------------------------------------------
------------------ ExperimentDef instance ------------------
------------------------------------------------------------

save :: BORL St Act -> IO ()
save borl = do
  res <- toSerialisableWith serializeSt id borl
  B.writeFile "saved.bin" (runPut $ put res)

load :: IO (Maybe (BORL St Act))
load = do
  bs <- liftIO $ B.readFile "saved.bin"
  case trace ("S.runGet") S.runGet S.get bs of
    Left err -> putStrLn ("Could not deserialize to Serializable BORL. Error: " ++ err) >> return Nothing
    Right ser -> trace ("ser") $ do
     borl <- buildBORLGrenade
     let (St sim _ _ _) = borl ^. s
     Just <$> fromSerialisableWith (deserializeSt (simRelease sim) (simDispatch sim) (simShipment sim) (simProcessingTimes $ simInternal sim)) id action (borl ^. B.actionFilter) netInp ser

instance ExperimentDef (BORL St Act) where
  type ExpM (BORL St Act) = IO
  type Serializable (BORL St Act) = BORLSerialisable StSerialisable Act
  serialisable = do
    res <- toSerialisableWith serializeSt id
    return res
  deserialisable ser =
    unsafePerformIO $ do
      borl <- buildBORLGrenade
      let (St sim _ _ _) = borl ^. s
      return $
        fromSerialisableWith
          (deserializeSt (simRelease sim) (simDispatch sim) (simShipment sim) (simProcessingTimes $ simInternal sim))
          id
          action
          (borl ^. B.actionFilter)
          netInp
          ser
  type InputValue (BORL St Act) = [Order]
  type InputState (BORL St Act) = ()
  -- ^ Generate some input values and possibly modify state. This function can be used to change the state. It is called
  -- before `runStep` and its output is used to call `runStep`.
  generateInput !_ !borl !_ !_ = do
    let !(St !_ !inc !_ !_) = borl ^. s
    return (inc, ())
  -- ^ Run a step of the environment and return new state and result.
  -- runStep :: (MonadIO m) => a -> InputValue a -> E.Period -> m ([StepResult], a)
  runStep !phase !borl !incOrds !_ = do
    !borl' <- stepM (set (s . nextIncomingOrders) incOrds borl)
    -- helpers
    when (borl ^. t `mod` 10000 == 0) $ liftIO $ prettyBORLMWithStInverse (Just $ mInverse borl) (setStPrettyPrintElems borl) >>= print
    let !simT = timeToDouble $ simCurrentTime $ borl' ^. s . simulation
    let !borlT = borl' ^. t
    -- demand
    let !demand = StepResult "demand" (Just simT) (fromIntegral $ length $ borl ^. s . nextIncomingOrders)
    -- cost related measures
    let (StatsOrderCost !earnOld !wipOld !boOld !fgiOld) = simStatsOrderCosts $ simStatistics (borl ^. s . simulation)
    let (StatsOrderCost !earn !wip !bo !fgi) = simStatsOrderCosts $ simStatistics (borl' ^. s . simulation)
    let !cEarn  = StepResult "EARN" (Just simT) (fromIntegral (earn - earnOld))
    let !cBoc   = StepResult "BOC" (Just simT) (boCosts costConfig * fromIntegral (bo - boOld))
    let !cWip   = StepResult "WIPC" (Just simT) (wipCosts costConfig * fromIntegral (wip - wipOld))
    let !cFgi   = StepResult "FGIC" (Just simT) (fgiCosts costConfig * fromIntegral (fgi - fgiOld))
    let !cSum   = StepResult "SUMC" (Just simT) (cBoc ^. resultYValue + cWip ^. resultYValue + cFgi ^. resultYValue)
    let !curOp  = StepResult "op" (Just simT) (fromIntegral $ length $ simOrdersOrderPool $ borl' ^. s . simulation)
    let !curWip = StepResult "wip" (Just simT) (fromIntegral $ wip - wipOld)
    let !curBo  = StepResult "bo" (Just simT) (fromIntegral $ bo - boOld)
    let !curFgi = StepResult "fgi" (Just simT) (fromIntegral $ fgi - fgiOld)
    -- time related measures
    let (StatsFlowTime !ftNrFloorAndFgi !(StatsOrderTime !sumTimeFloorAndFgi !stdDevFloorAndFgi !_) !mTardFloorAndFgi) = simStatsShopFloorAndFgi $ simStatistics (borl' ^. s . simulation)
    let !tFtMeanFloorAndFgi     = StepResult "FTMeanFloorAndFgi" (Just simT) (fromRational sumTimeFloorAndFgi / fromIntegral ftNrFloorAndFgi)
    let !tFtStdDevFloorAndFgi   = StepResult "FTStdDevFloorAndFgi" (Just simT) (maybe 0 fromRational $ getWelfordStdDev stdDevFloorAndFgi)
    let !tTardPctFloorAndFgi    = StepResult "TARDPctFloorAndFgi" (Just simT) (maybe 0 (\(StatsOrderTard nrTard sumTard stdDevTard) -> fromIntegral nrTard / fromIntegral ftNrFloorAndFgi) mTardFloorAndFgi)
    let !tTardMeanFloorAndFgi   = StepResult "TARDMeanFloorAndFGI" (Just simT) (maybe 0 (\(StatsOrderTard nrTard sumTard stdDevTard) -> fromRational sumTard / fromIntegral nrTard) mTardFloorAndFgi)
    let !tTardStdDevFloorAndFgi = StepResult "TARDStdDevFloorAndFGI" (Just simT) (maybe 0 (\(StatsOrderTard nrTard sumTard stdDevTard) -> maybe 0 fromRational $ getWelfordStdDev stdDevTard) mTardFloorAndFgi)
    let !(StatsFlowTime !ftNrFloor !(StatsOrderTime !sumTimeFloor !stdDevFloor !_) !mTardFloor) = simStatsShopFloor $ simStatistics (borl' ^. s . simulation)
    let !tFtMeanFloor     = StepResult "FTMeanFloor" (Just simT) (fromRational sumTimeFloor / fromIntegral ftNrFloor)
    let !tFtStdDevFloor   = StepResult "FTStdDevFloor" (Just simT) (maybe 0 fromRational $ getWelfordStdDev stdDevFloor)
    let !tTardPctFloor    = StepResult "TARDPctFloor" (Just simT) (maybe 0 (\(StatsOrderTard nrTard sumTard stdDevTard) -> fromIntegral nrTard / fromIntegral ftNrFloor) mTardFloor)
    let !tTardMeanFloor   = StepResult "TARDMeanFloor" (Just simT) (maybe 0 (\(StatsOrderTard nrTard sumTard stdDevTard) -> fromRational sumTard / fromIntegral nrTard) mTardFloor)
    let !tTardStdDevFloor = StepResult "TARDStdDevFloor" (Just simT) (maybe 0 (\(StatsOrderTard nrTard sumTard stdDevTard) -> maybe 0 fromRational $ getWelfordStdDev stdDevTard) mTardFloor)
    -- BORL' related measures
    let valS !l = realToFrac (borl' ^?! l)
        valV !l = realToFrac $ head $ fromValue (borl' ^?! l)
        valL !l = realToFrac $ V.head $ borl' ^?! l
    let !avgRew    = StepResult "AvgReward" (Just $ fromIntegral borlT) (valL $ proxies . rho . proxyScalar)
        !expAvgRew = StepResult "ExpAvgReward" (Just $ fromIntegral borlT) (valS $ expSmoothedReward)
        !avgRewMin = StepResult "MinAvgReward" (Just $ fromIntegral borlT) (valL $ proxies . rhoMinimum . proxyScalar)
        !pltP1     = StepResult "PLT P1" (Just $ fromIntegral borlT) (timeToDouble $ M.findWithDefault 0 (Product 1) (borl' ^. s . plannedLeadTimes))
        !pltP2     = StepResult "PLT P2" (Just $ fromIntegral borlT) (timeToDouble $ M.findWithDefault 0 (Product 2) (borl' ^. s . plannedLeadTimes))
        !psiRho    = StepResult "PsiRho" (Just $ fromIntegral borlT) (valV $ psis . _1)
        !psiV      = StepResult "PsiV" (Just $ fromIntegral borlT) (valV $ psis . _2)
        !psiW      = StepResult "PsiW" (Just $ fromIntegral borlT) (valV $ psis . _3)
        !vAvg      = StepResult "VAvg" (Just $ fromIntegral borlT) (avg $ V.toList $ borl' ^. lastRewards)
        !reward    = StepResult "Reward" (Just $ fromIntegral borlT) (realToFrac $ headWithDefault 0 $ borl' ^. lastRewards)
        avg !xs    = realToFrac $ sum xs / fromIntegral (length xs)
        headWithDefault d vec    = fromMaybe d $ vec V.!? 0
    return $! force $
      if phase /= EvaluationPhase
      then ([
        --   avgRew
        -- , avgRewMin
        -- , expAvgRew
        -- -- , vAvg
        -- , reward

        ], borl')

      else
      ( [-- cost related measures
          cSum
        , cEarn
        , cBoc
        , cWip
        , cFgi
             -- floor
        , curOp
        , curWip
        , curBo
        , curFgi
        , demand
             -- time related measures
        , tFtMeanFloorAndFgi
        , tFtStdDevFloorAndFgi
        , tTardPctFloorAndFgi
        , tTardMeanFloorAndFgi
        , tTardStdDevFloorAndFgi
        , tFtMeanFloor
        , tFtStdDevFloor
        , tTardPctFloor
        , tTardMeanFloor
        , tTardStdDevFloor
             -- BORL related measures
        , avgRew
        -- , avgRewMin
        -- , expAvgRew
        -- , pltP1
        -- , pltP2
        -- , psiRho
        -- , psiV
        -- , psiW
        -- , vAvg
        -- , reward
        ]
      , borl')
  -- ^ Provides the parameter setting.
  -- parameters :: a -> [ParameterSetup a]
  parameters borl =
    [ ParameterSetup
        "Algorithm"
        (set algorithm)
        (view algorithm)
        (Just $ return .
         const
           [
             AlgDQNAvgRewAdjusted 0.8 1.0 ByStateValues
           , AlgDQNAvgRewAdjusted 0.8 0.99 ByStateValues
           , AlgDQN 0.99 Exact
           ])
        Nothing
        Nothing
        Nothing
    , ParameterSetup
        "RewardType"
        (set (s . rewardFunctionOrders))
        (view (s . rewardFunctionOrders))
        (Just $ return .
         const
           [
             -- RewardPeriodEndSimple (ConfigRewardCosts (Just 750))
             RewardPeriodEndSimple (ConfigRewardCosts (Just 200))
           ])
        Nothing
        Nothing
        Nothing
    , ParameterSetup
        "ReleaseAlgorithm"
        (\r -> over (s . simulation) (\sim -> sim {simRelease = r}))
        (simRelease . view (s . simulation))
        (Just . const . return $
           ([ mkReleasePLT initialPLTS
            -- , releaseImmediate
            -- , releaseBIL (M.fromList [(Product 1, 6), (Product 2, 6)])
            -- , releaseBIL (M.fromList [(Product 1, 5), (Product 2, 5)])
            -- , releaseBIL (M.fromList [(Product 1, 4), (Product 2, 4)])
            -- , releaseBIL (M.fromList [(Product 1, 3), (Product 2, 3)])
            -- , releaseBIL (M.fromList [(Product 1, 2), (Product 2, 2)])
            -- , releaseBIL (M.fromList [(Product 1, 1), (Product 2, 1)])
           ] ++
            map (\lts -> releaseBIL (M.fromList (map (\pt -> (pt, lts)) productTypes))) [1..4]
        ))
        Nothing
        (Just (\x -> uniqueReleaseName x /= pltReleaseName)) -- drop preparation phase for all release algorithms but the BORL releaser
        (Just
           (\x ->
              if uniqueReleaseName x == pltReleaseName
                then FullFactory
                else SingleInstance -- only evaluate once if ImRe or BIL
            ))
    ] ++
    -- [ ParameterSetup
    --   "Xi (at period 0)"
    --   (set (B.parameters . xi))
    --   (^. B.parameters . xi)
    --   (Just $ return . const [5e-3])
    --   Nothing Nothing Nothing
    -- ] ++
    -- [ ParameterSetup
    --   "Zeta (at period 0)"
    --   (set (B.parameters . zeta))
    --   (^. B.parameters . zeta)
    --   (Just $ return . const [0.10])
    --   Nothing Nothing Nothing
    -- ] ++
    [ ParameterSetup
      "Alpha (at period 0)"
      (set (B.parameters . alpha))
      (^. B.parameters . alpha)
      (Just $ return . const [0.01])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "AlphaRhoMin (at period 0)"
      (set (B.parameters . alphaRhoMin))
      (^. B.parameters . alphaRhoMin)
      (Just $ return . const [2e-5])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "Delta (at period 0)"
      (set (B.parameters . delta))
      (^. B.parameters . delta)
      (Just $ return . const [0.005])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "Gamma (at period 0)"
      (set (B.parameters . gamma))
      (^. B.parameters . gamma)
      (Just $ return . const [0.01])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "Epsilon (at period 0)"
      (set (B.parameters . epsilon))
      (^. B.parameters . epsilon)
      (Just $ return . const [fromList [0.25, 0.25]])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "Exploration (at period 0)"
      (set (B.parameters . exploration))
      (^. B.parameters . exploration)
      (Just $ return . const [1.0])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "Learn Random Above (faster converging rho)"
      (set (B.parameters . learnRandomAbove))
      (^. B.parameters . learnRandomAbove)
      (Just $ return . const [0.5])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "Decay Alpha"
      (set (B.decaySetting . alpha))
      (^. B.decaySetting . alpha)
      (Just $ return . const [ExponentialDecay (Just 5e-5) 0.25 50000])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "Decay AlphaRhoMin"
      (set (B.decaySetting . alphaRhoMin))
      (^. B.decaySetting . alphaRhoMin)
      (Just $ return . const [ExponentialDecay (Just 2e-5) 0.25 50000])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "Decay Delta"
      (set (B.decaySetting . delta))
      (^. B.decaySetting . delta)
      (Just $ return . const [ExponentialDecay (Just 5e-4) 0.25 100000])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "Decay Gamma"
      (set (B.decaySetting . gamma))
      (^. B.decaySetting . gamma)
      (Just $ return . const [ExponentialDecay (Just 1e-3) 0.25 100000])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "Decay Epsilon"
      (set (B.decaySetting . epsilon))
      (^. B.decaySetting . epsilon)
      (Just $ return . const [fromList [NoDecay]])
      Nothing Nothing Nothing
    ] ++
    [ ParameterSetup
      "Decay Exploration"
      (set (B.decaySetting . exploration))
      (^. B.decaySetting . exploration)
      (Just $ return . const [ExponentialDecay (Just 0.005) 0.25 100000])
      Nothing Nothing Nothing
    ] ++
    -- Cannot be changed here!!!
    -- [ ParameterSetup
    --   "Replay Memory Size"  Cannot be changed here as initialised!!!
    --   (setAllProxies (proxyNNConfig . replayMemoryMaxSize))
    --   (^?! proxies . v . proxyNNConfig . replayMemoryMaxSize)
    --   (Just $ return . const [10000])
    --   Nothing
    --   Nothing
    --   Nothing
    -- | isNN
    -- ] ++
    -- [ ParameterSetup
    --   "Replay Memory Strategy" Cannot be changed here as initialised!!!!!!
    --   (setAllProxies (proxyNNConfig . replayMemoryStrategy))
    --   (^?! proxies . v . proxyNNConfig . replayMemoryStrategy)
    --   (Just $ return . const [ReplayMemoryPerAction])
    --   Nothing
    --   Nothing
    --   Nothing
    -- | isNN
    -- ] ++
    [ ParameterSetup
      "Training Batch Size"
      (setAllProxies (proxyNNConfig . trainBatchSize))
      (^?! proxies . v . proxyNNConfig . trainBatchSize)
      (Just $ return . const [4])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "Training Iterations"
      (setAllProxies (proxyNNConfig . trainingIterations))
      (^?! proxies . v . proxyNNConfig . trainingIterations)
      (Just $ return . const [1])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "ANN (Grenade) Learning Rate"
      (setAllProxies (proxyNNConfig . grenadeLearningParams))
      (^?! proxies . v . proxyNNConfig . grenadeLearningParams)
      (Just $ return . const [OptAdam 0.005 0.9 0.999 1e-8 1e-3 -- , OptAdam 0.0001 0.9 0.999 1e-8 1e-3
                             ])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "ANN Learning Rate Decay"
      (setAllProxies (proxyNNConfig . learningParamsDecay))
      (^?! proxies . v . proxyNNConfig . learningParamsDecay)
      (Just $ return . const [ ExponentialDecay (Just 5e-6) (configDecayRate decay) (configDecaySteps decay)])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "Grenade Smooth Target Update"
      (setAllProxies (proxyNNConfig . grenadeSmoothTargetUpdate))
      (^?! proxies . v . proxyNNConfig . grenadeSmoothTargetUpdate)
      (Just $ return . const [0.01])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "Grenade Smooth Target Update Period"
      (setAllProxies (proxyNNConfig . grenadeSmoothTargetUpdatePeriod))
      (^?! proxies . v . proxyNNConfig . grenadeSmoothTargetUpdatePeriod)
      (Just $ return . const [100])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "ScaleParameters"
      (setAllProxies (proxyNNConfig . scaleParameters))
      (^?! proxies . v . proxyNNConfig . scaleParameters)
      (Just $ return . const [ScalingNetOutParameters (-800) 800 (-300) 300 (-400) 800 (-400) 800])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "Scaling Out Setup"
      (setAllProxies (proxyNNConfig . scaleOutputAlgorithm))
      (^?! proxies . v . proxyNNConfig . scaleOutputAlgorithm)
      (Just $ return . const [ScaleMinMax])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "Grenade Dropout Layers Flip-Active Periods"
      (setAllProxies (proxyNNConfig . grenadeDropoutFlipActivePeriod))
      (^?! proxies . v . proxyNNConfig . grenadeDropoutFlipActivePeriod)
      (Just $ return . const [10^5])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "Grenade Dropout Only-Active After Period"
      (setAllProxies (proxyNNConfig . grenadeDropoutFlipActivePeriod))
      (^?! proxies . v . proxyNNConfig . grenadeDropoutFlipActivePeriod)
      (Just $ return . const [0])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "Gradient clipping"
      (setAllProxies (proxyNNConfig . clipGradients))
      (^?! proxies . v . proxyNNConfig . clipGradients)
      (Just $ return . const [NoClipping])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "NStep"
      (set (settings . nStep))
      (^. settings . nStep)
      (Just $ return . const [5])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "Exploration Strategy"
      (set (settings . explorationStrategy))
      (^. settings . explorationStrategy)
      (Just $ return . const [EpsilonGreedy])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "Workers Min Exploration"
      (set (settings . workersMinExploration))
      (^. settings . workersMinExploration)
      (Just $ return . const [take 8 $ 0.01 : 0.02 : 0.03 : 0.04 : [0.05, 0.10 .. 1.0]])
      Nothing
      Nothing
      Nothing
    | isNN
    ] ++
    [ ParameterSetup
      "Independent Agents Share Rhos"
      (set (settings . independentAgentsSharedRho))
      (^. settings . independentAgentsSharedRho)
      (Just $ return . const [True])
      Nothing
      Nothing
      Nothing
    ] ++
    [ ParameterSetup
      "Overestimate Rho"
      (set (settings . overEstimateRho))
      (^. settings . overEstimateRho)
      (Just $ return . const [False])
      Nothing
      Nothing
      Nothing
    ] ++
    [ ParameterSetup
      "Main Agent Selects Greedy Actions"
      (set (settings . mainAgentSelectsGreedyActions))
      (^. settings . mainAgentSelectsGreedyActions)
      (Just $ return . const [False])
      Nothing
      Nothing
      Nothing
    ]
    where
      isNN = isNeuralNetwork (borl ^. proxies . v)
  -- HOOKS
  -- beforePreparationHook _ _ g borl =
  --   liftIO $ do
  --     let dir = "results/" <> T.unpack (T.replace " " "_" experimentName) <> "/data/"
  --     createDirectoryIfMissing True dir
  --     writeFile (dir ++ "plot.sh") gnuplot
  --     mapMOf (s . simulation) (setSimulationRandomGen g) borl
  beforeWarmUpHook _ _ _ g borl = liftIO $ mapMOf (s . simulation) (fmap resetStatistics . setSimulationRandomGen g) $ set (B.parameters . exploration) 0.00 $ set (B.settings . disableAllLearning) True borl
  beforeEvaluationHook _ _ _ g borl = liftIO $

    mapMOf (s . simulation) (fmap resetStatistics . setSimulationRandomGen g) $ set (B.parameters . exploration) 0.00 $ set (B.settings . disableAllLearning) True  borl
  -- afterPreparationHook _ expNr repetNr = liftIO $ copyFiles "prep_" expNr repetNr Nothing
  -- afterWarmUpHook _ expNr repetNr repliNr = liftIO $ copyFiles "warmup_" expNr repetNr (Just repliNr)
  -- afterEvaluationHook _ expNr repetNr repliNr = liftIO $ copyFiles "eval_" expNr repetNr (Just repliNr)

setStPrettyPrintElems :: BORL St Act -> BORL St Act
setStPrettyPrintElems borl = setAllProxies (proxyNNConfig . prettyPrintElems) [netInp $ borl ^. s] borl

expSetting :: BORL St Act -> ExperimentSetting
expSetting borl =
  ExperimentSetting
    { _experimentBaseName = experimentName
    , _experimentInfoParameters = [actBounds, pltBounds, csts, dem, ftExtr, rout, procT, isNN] ++ concat [[repMemSize, repMemStrat, nnSetup] | isNNFlag]
    , _experimentRepetitions = 1
    , _preparationSteps = 300000
    , _evaluationWarmUpSteps = 1000
    , _evaluationSteps = 7000
    , _evaluationReplications = 30
    , _evaluationMaxStepsBetweenSaves = Just 6000
    }
  where
    isNNFlag = isNeuralNetwork (borl ^. proxies . v)
    isNN = ExperimentInfoParameter "Is Neural Network" isNNFlag
    -- dec = ExperimentInfoParameter "Decay" (configDecayName decay) -- decay is not used in experiment
    actBounds = ExperimentInfoParameter "Action Bounds" (configActLower actionConfig, configActUpper actionConfig)
    pltBounds = ExperimentInfoParameter "Action Filter (Min/Max PLT)" (configActFilterMin actionFilterConfig, configActFilterMax actionFilterConfig)
    csts = ExperimentInfoParameter "Costs" costConfig
    dem = ExperimentInfoParameter "Demand" (configDemandName demand)
    ftExtr = ExperimentInfoParameter "Feature Extractor (State Representation)" (configFeatureExtractorName $ featureExtractor True)
    rout = ExperimentInfoParameter "Routing (Simulation Setup)" (configRoutingName routing)
    procT = ExperimentInfoParameter "Processing Time (Simulation Setup)" (configProcTimesName procTimes)
    repMemSize = ExperimentInfoParameter "Replay Memory Size" (borl ^?! proxies . v . proxyNNConfig . replayMemoryMaxSize)
    repMemStrat = ExperimentInfoParameter "Replay Memory Strategy" (borl ^?! proxies . v . proxyNNConfig . replayMemoryStrategy)
    nnSetup = ExperimentInfoParameter "ANN Architecture" (netSpec $ borl ^?! proxies . v)
    netSpec px = case px of
      Grenade t _ _ _ _ _ -> T.replace "Spec" "" $ T.replace ",1,1)" ")" $ T.replace "SpecRelu " "Relu" $ T.replace "SpecFullyConnected" "FC" $ tshow (networkToSpecification t)
      CombinedProxy s _ _ -> netSpec s
      _ -> ("" :: T.Text)

max_stars_repo_name: schnecki/borl-releaser

max_issues_repo_name: schnecki/borl-releaser

max_forks_repo_name: schnecki/borl-releaser


text:
module Quantum where

import Data.AdditiveGroup
import Data.Complex hiding ((:+))
import qualified Data.Complex as C (Complex ((:+)))
import Data.Foldable (foldl')
import Data.VectorSpace hiding (magnitude)
import Lib
import Matrix
import Vector

-- base type for qubits
type CDouble = Complex Double

type CVVec n = VVec n CDouble

type Qubit = CVVec Two

type NQubits n = CVVec (Exp Two n)

-- for ease of use
sqtwo :: Floating a => a
sqtwo = 1 / sqrt 2

i :: Num a => Complex a
i = 0 C.:+ 1

-- define some basic qubits/vectors, and qubit collections
basisVec :: forall n a. (Num a, KnownNat n) => Fin n -> VVec n a
basisVec f = generateMat (\x _ -> fromIntegral $ fromEnum $ x == f)

zero' :: forall n a. (Num a, KnownNat n) => VVec n a
zero' = basisVec FZero

zero :: Qubit
zero = zero'

one' :: forall n a m. (Num a, KnownNat n, n ~ 'Succ m) => VVec n a
one' = basisVec (FSucc FZero)

one :: Qubit
one = one'

plus :: Qubit
plus = sqtwo *^ zero ^+^ sqtwo *^ one

minus :: Qubit
minus = sqtwo *^ zero ^-^ sqtwo *^ one

qubits :: [Qubit]
qubits = [zero, one, plus, minus]

makeTwoQubits :: [Qubit] -> [NQubits Two]
makeTwoQubits qs = (.*.) <$> qs <*> qs

twoQubits :: [NQubits Two]
twoQubits = makeTwoQubits [zero, one]

consQubit :: CDouble -> CDouble -> Qubit
consQubit a b = numMatFromList [[a], [b]]

consVVec :: forall n a. (KnownNat n, Num a) => [a] -> VVec n a
consVVec ns = numMatFromList @n @One (fmap (: []) ns)

-- define some transformation matrices
pauliX :: Num a => Matrix Two Two a
pauliX = numMatFromList [[0, 1], [1, 0]]

pauliY :: RealFloat a => Matrix Two Two (Complex a)
pauliY = numMatFromList [[0, - i], [i, 0]]

pauliZ :: Num a => Matrix Two Two a
pauliZ = numMatFromList [[1, 0], [0, -1]]

hadamard :: Floating a => Matrix Two Two a
hadamard = (* sqtwo) <$> numMatFromList [[1, 1], [1, -1]]

-- a rotation matrix. produces a matrix which rotates by n radians
rotation :: Floating a => a -> Matrix Two Two a
rotation n = numMatFromList [[c, - s], [s, c]]
  where
    c = cos n
    s = sin n

cnot :: Num a => Matrix Four Four a
cnot = numMatFromList [[1, 0, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]

cnot' :: Num a => Matrix Four Four a
cnot' = numMatFromList [[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]]

iden2 :: Num a => Matrix Two Two a
iden2 = identity

swap :: Num a => Matrix Four Four a
swap = numMatFromList [[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]]

-- define the tensor product
tensorProd ::
  (Num a) => Matrix n m a -> Matrix i j a -> Matrix (Mul n i) (Mul m j) a
tensorProd n m = expandNested $ fmap (\a -> fmap (a *) m) n

-- ease of use function for the tensor product
(.*.) :: (Num a) => Matrix n m a -> Matrix i j a -> Matrix (Mul n i) (Mul m j) a
(.*.) = tensorProd

infixr 8 .*.

-- apply the tensor product a number of times.
-- to use this function, use `tensorPower @(type number) matrix`
-- for example, to tensor product matrix `a` gate three times, do `tensorPower @Three a`
-- If used in relation to other matrices, tensorPower can infer
--  the `i` value, which is pretty neat
-- Only works on non-vector matrices - see tensorPowerVVec for column vec variant
tensorPower ::
  forall i a n m.
  (Num a, KnownNat i, i ~ GetExp ('Succ n) i, i ~ GetExp ('Succ m) i) =>
  Matrix ('Succ n) ('Succ m) a ->
  Matrix (Exp ('Succ n) i) (Exp ('Succ m) i) a
tensorPower m = tensorPower' @i m natSing

-- apply the tensor product a number of times on a column vector
-- very similar to tensorPower, but only works on column vectors
tensorPowerVVec ::
  forall i a n.
  (Num a, KnownNat i, i ~ GetExp ('Succ n) i) =>
  Matrix ('Succ n) 'One a ->
  Matrix (Exp ('Succ n) i) (Exp 'One i) a
tensorPowerVVec m = tensorPower' @i m natSing

-- helper function for tensorPower
tensorPower' ::
  forall i n m a.
  Num a =>
  Matrix n m a ->
  NatS i ->
  Matrix (Exp n i) (Exp m i) a
tensorPower' m OneS = m
tensorPower' m (SuccS s) = m .*. tensorPower' m s

-- transpose and get the conjugate of the given qubit
-- effectively, find the bra
conjTrans :: CVVec n -> HVec n CDouble
conjTrans v = conjugate <$> Matrix.transpose v

-- using conjTrans, find the inner product of two qubits
innerProduct :: CVVec n -> CVVec n -> CDouble
innerProduct v v' = getVal $ conjTrans v *.* v'

-- generate all basis vectors for a given dimension
compBasis' :: forall n a. (Num a, KnownNat n) => Vector n (VVec n a)
compBasis' = generateVec basisVec

-- define the computational basis
compBasis :: Vector Two Qubit
compBasis = compBasis'

-- apply a transformation matrix to the computational basis
transformBasis ::
  forall n a. (Num a, KnownNat n) => Matrix n n a -> Vector n (VVec n a)
transformBasis m = fmap (m *.*) compBasis'

-- measure in a given basis, returning probabilities for each
-- basis vector
measureIn :: CVVec m -> Vector n (CVVec m) -> Vector n Double
measureIn v = fmap ((** 2) . magnitude . Quantum.innerProduct v)

-- measure in the computational basis, returning probabilities
-- for each basis vector
measure :: Qubit -> Vector Two Double
measure v = measureIn v compBasis

-- given a list of vector-probability pairs, return Just the density matrix
--  if the probabilities sum up to 1; else return Nothing
getDensityMatrix ::
  (KnownNat m) => [(CVVec m, CDouble)] -> Maybe (Matrix m m CDouble)
getDensityMatrix qs
  | sum (map snd qs) == 1 =
    Just $
      foldr (^+^) zeroed [p *^ f tq | (tq, p) <- qs]
  | otherwise = Nothing
  where
    zeroed = generateMat (\_ _ -> 0)
    f m = m *.* conjTrans m

-- given a density matrix (as seen in previous function),
--  and a complex number vector, return the probability of
--  seeing that vector
calcProbability :: Matrix m m CDouble -> CVVec m -> CDouble
calcProbability densityMatrix base =
  getVal $ conjTrans base *.* densityMatrix *.* base

chi_f :: CVVec m -> CVVec m -> CDouble
chi_f s x = (-1) ** Quantum.innerProduct s x

-- vecToFunc :: (KnownNat m) => CVVec m -> Matrix m m CDouble
-- vecToFunc v =
--   generateMat (\f f' -> if f == f' then getAtMatrix f FZero v else 0)

-- fhat :: (KnownNat n) => (CVVec n -> CDouble) -> CVVec n -> CDouble
-- fhat f s =
--   (foldl' (\b a -> b + (apply a)) 0 compBasis') / (2 ** (fst $ Matrix.size s))
--  where
--   chi_s = chi s
--   apply x = (f x) * chi_s x

createInput :: Num a => Vector n (VVec m a) -> VVec (Exp m n) a
createInput (VecSing v) = v
createInput (v :+ vs) = v .*. createInput vs

fourierBasis :: (KnownNat n) => Vector n (CVVec n)
fourierBasis = fmap chi compBasis'

chi :: (KnownNat n) => CVVec n -> CVVec n
chi s = normalise $ toVVec $ fmap (chi_f s) compBasis'

fhat_f :: (KnownNat n) => CVVec n -> CVVec n -> CDouble
fhat_f = Quantum.innerProduct

fhat :: (KnownNat n) => CVVec n -> CVVec n
fhat f = toVVec $ fmap (fhat_f f) fourierBasis

normalise :: (KnownNat n) => CVVec n -> CVVec n
normalise v = v ^/ sqrt (sum $ fmap (abs) v)

groverDiffusion ::
  forall i a.
  (KnownNat (Exp Two i), KnownNat i, i ~ GetExp Two i, Floating a) =>
  Matrix (Exp Two i) (Exp Two i) a
groverDiffusion =
  tensorPower hadamard *.* generateMat f *.* tensorPower hadamard
  where
    f a b
      | a == b && a == FZero = 1
      | a == b = -1
      | otherwise = 0

-- fhat :: (KnownNat n) => CVVec n -> CVVec n -> CDouble
-- fhat f s

ext: hs


text:
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE BangPatterns,
             ScopedTypeVariables,
             RecordWildCards,
             FlexibleContexts,
             TypeFamilies #-}
import Data.IDX (IDXData,  decodeIDXFile, idxIntContent )
{------------------------}
import Numeric.LinearAlgebra as NL (cmap, (#>), vjoin, Numeric,  Vector, fromList )
{------------------------}
import Data.Vector.Unboxed as DTU ( toList )
import qualified Data.Vector.Storable as DT
import Data.List (sort)
import Prelude hiding (readFile)
import Neuro
    ( ActivationFunction, Network(..), loadNetwork )

kol_vo :: Int
kol_vo = 1000

main :: IO ()
main = do
        res <- loadNetwork "smartNet5.nn"
        Just idxTrain  <- decodeIDXFile "test/t10k-images.idx3-ubyte" -- image
        Just idxResult <- decodeIDXFile "test/t10k-labels.idx1-ubyte" -- result
        --writeFile "resNeuro.txt" (concat $ take (28*kol_vo) $ drowImage idxTrain)
        --writeFile "resLabel.txt" (show $ take kol_vo $ DTU.toList $  idxIntContent idxResult)
        putStrLn "Точность алгоритма :"
        print $ (\x -> (fromIntegral x) / fromIntegral kol_vo) $ length $ filter (\((_,x),y) -> x == y) $ flip zip (filterCorrect idxResult) $  printRes kol_vo res idxTrain $ DTU.toList $ idxIntContent idxResult
        where filterCorrect idxResult = reverse $ take kol_vo $ DTU.toList $ idxIntContent idxResult

printRes :: Int -> Network Double -> IDXData -> t -> [(Double,Int)]
printRes 0 _   _           _      = []
printRes x res idxTrain idxResult = last (sort $ flip zip [0..] $ softmax $ DT.toList $ Main.output res tanh (last $ take x $ image idxTrain)) : printRes (x-1) res idxTrain idxResult

output :: (Floating (Vector a), Numeric a, DT.Storable a, Num (Vector a)) => Network a -> ActivationFunction a -> Vector a -> Vector a
output (Network{..}) act input = foldl f (vjoin [input, 1]) matrices
  where f (!inp) m = cmap act $ m #> inp

softmax :: [Double] -> [Double]
softmax xs = let xs' = exp <$> xs
                 s   = sum xs'
             in map (/ s) xs'

drowImage :: IDXData -> [String]
drowImage = _0_255to0or1 . DTU.toList . idxIntContent

_0_255to0or1 :: [Int] -> [String]
_0_255to0or1 [] = []
_0_255to0or1 xs = (show (concat $ fmap show $ convert $ take 28 xs) ++ " \n") : _0_255to0or1 (drop 28 xs)

convert [] = []
convert (x:xs)
              | x > 240 = 1 : convert xs
              | otherwise = 0 : convert xs

image :: IDXData -> [Vector Double]
image  = matrix2x2 . _0_1to0_255 . DTU.toList . idxIntContent
result :: IDXData -> [Vector Double]
result = unitar . DTU.toList . idxIntContent

unitar :: [Int] -> [NL.Vector Double]
unitar []     = []
unitar (0:xs) = NL.fromList [1,-1,-1,-1,-1,-1,-1,-1,-1,-1] : unitar xs
unitar (1:xs) = NL.fromList [-1,1,-1,-1,-1,-1,-1,-1,-1,-1] : unitar xs
unitar (2:xs) = NL.fromList [-1,-1,1,-1,-1,-1,-1,-1,-1,-1] : unitar xs
unitar (3:xs) = NL.fromList [-1,-1,-1,1,-1,-1,-1,-1,-1,-1] : unitar xs
unitar (4:xs) = NL.fromList [-1,-1,-1,-1,1,-1,-1,-1,-1,-1] : unitar xs
unitar (5:xs) = NL.fromList [-1,-1,-1,-1,-1,1,-1,-1,-1,-1] : unitar xs
unitar (6:xs) = NL.fromList [-1,-1,-1,-1,-1,-1,1,-1,-1,-1] : unitar xs
unitar (7:xs) = NL.fromList [-1,-1,-1,-1,-1,-1,-1,1,-1,-1] : unitar xs
unitar (8:xs) = NL.fromList [-1,-1,-1,-1,-1,-1,-1,-1,1,-1] : unitar xs
unitar (9:xs) = NL.fromList [-1,-1,-1,-1,-1,-1,-1,-1,-1,1] : unitar xs

_0_1to0_255 :: [Int] -> [Double]
_0_1to0_255 [] = []
_0_1to0_255 (x:xs) = (fromIntegral x / 255) : _0_1to0_255 xs

matrix2x2 :: [Double] -> [NL.Vector Double]
matrix2x2 [] = []
matrix2x2 xs = NL.fromList (Prelude.take 784 xs) : matrix2x2 (Prelude.drop 784 xs)


text:
module STCR2S1PointSet where

import           Control.Monad           as M
import           Control.Monad.Parallel  as MP
import           Data.Array.Repa         as R
import           Data.Binary             (decodeFile)
import           Data.Complex
import           Data.List               as L
import           Data.Vector.Unboxed     as VU
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.MonteCarlo
import           Image.IO
import           STC.CompletionFieldR2S1
import STC.Shape
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Types
import           Utils.Array
import FokkerPlanck.Histogram
import FokkerPlanck.GreensFunction

main = do
  args@(numPointStr:numOrientationStr:sigmaStr:taoStr:numTrailStr:maxTrailStr:initOriStr:initSpeedStr:writeFlagStr:numIterationStr:thresholdStr:histFilePath:rStr:deltaStr:stdStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      sigma = read sigmaStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      initOri = read initOriStr :: Double
      initSpeed = read initSpeedStr :: Double
      writeFlag = read writeFlagStr :: Bool
      numIteration = read numIterationStr :: Int
      threshold = read thresholdStr :: Double
      r = read rStr :: Double
      delta = read deltaStr :: Double
      std = read stdStr :: Double
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCR2S1PointSet"
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  -- arrG' <-
  --   if flag
  --     then getNormalizedHistogramArr <$> decodeFile histFilePath
  --     else do
  --       putStrLn "Couldn't find a Green's function data. Start simulation..."
  --       solveMonteCarloR2S1
  --         numThread
  --         numTrail
  --         maxTrail
  --         numPoint
  --         numPoint
  --         numOrientation
  --         sigma
  --         tao
  --         ((sqrt 2) * fromIntegral numPoint)
  --         initSpeed
  --         histFilePath
  --         -- (0, 0, 0, 0, initOri / 180 * pi, initSpeed)
  arrG' <-
    sampleR2S1
      numPoint
      numPoint
      delta
      (2 * pi / fromIntegral numOrientation)
      initSpeed
      sigma
      tao
      [0 .. numOrientation - 1]
  let -- arrG' =
      --   sampleR2S1
      --     numPoint
      --     numPoint
      --     delta
      --     (2 * pi / fromIntegral numOrientation)
      --     initSpeed
      --     sigma
      --     tao
      --     [0 .. numOrientation - 1]
      -- arrG =
      --   computeS . R.traverse arrG' id $ \f idx@(Z :. k :. i :. j) ->
      --     let r' =
      --           sqrt . fromIntegral $
      --           (i - div numPoint 2) ^ 2 + (j - div numPoint 2) ^ 2
      --     in -- if r' > 49 -- || r' <= 2
      --        --   then 0
      --        --   else
      --       f idx
      arrG = computeS . rotate3DR $ arrG'
      -- r = 16
  let numTheta = 8
      deltaTheta = (2 * pi) / numTheta
      xs -- [R2S1RPPoint (-i, 0, 0, 1) | i <- [-24,-18 .. 24]]
       =
        ((L.map
            (\(i, j) -> R2S1RPPoint (round i, round j, 0, 1))
            [ (r * cos (k * deltaTheta) + 0, r * sin (k * deltaTheta) + 0)
            | k <- [0 .. numTheta - 1]
            ]) -- L.++
         -- [R2S1RPPoint (-3, -1, 0, 1)]
         )
      bias = computeBias numPoint numPoint numOrientation xs
      -- initialEigenVec =
      --   computeInitialEigenVec numPoint numPoint numOrientation xs
      circle = Points (0, 0) 2 Circle {circleNum = 8, circleRadius = round r}
      circlePoints = makeShape2D $ circle
      ys = getShape2DIndexListGaussian 5 2 $ circlePoints
      initialEigenVec =
        computeInitialEigenVecGaussian numPoint numPoint numOrientation ys
  print xs
  print circlePoints
  -- powerMethod
  --   emptyPlan
  --   folderPath
  --   arrG
  --   numIteration
  --   writeFlag
  --   ""
  --   threshold
  --   bias
  --   initialEigenVec
  computeContourR2S1Analytic
    folderPath
    arrG
    numOrientation
    numPoint
    numPoint
    sigma
    tao
    initSpeed
    std
    delta
    bias
    circlePoints
  -- computeContourR2S1Tangent
  --   folderPath
  --   arrG
  --   numOrientation
  --   numPoint
  --   numPoint
  --   delta
  --   circlePoints

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE FlexibleContexts #-}
-- | Pearson's chi squared test.
module Statistics.Test.ChiSquared (
    chi2test
    -- * Data types
  , TestType(..)
  , TestResult(..)
  ) where

import Prelude hiding (sum)
import Statistics.Distribution
import Statistics.Distribution.ChiSquared
import Statistics.Function (square)
import Statistics.Sample.Internal (sum)
import Statistics.Test.Types
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U

-- | Generic form of Pearson chi squared tests for binned data. Data
--   sample is supplied in form of tuples (observed quantity,
--   expected number of events). Both must be positive.
chi2test :: (G.Vector v (Int,Double), G.Vector v Double)
         => Double              -- ^ p-value
         -> Int                 -- ^ Number of additional degrees of
                                --   freedom. One degree of freedom
                                --   is due to the fact that the are
                                --   N observation in total and
                                --   accounted for automatically.
         -> v (Int,Double)      -- ^ Observation and expectation.
         -> TestResult
chi2test p ndf vec
  | ndf < 0        = error $ "Statistics.Test.ChiSquare.chi2test: negative NDF " ++ show ndf
  | n   < 0        = error $ "Statistics.Test.ChiSquare.chi2test: too short data sample"
  | p > 0 && p < 1 = significant $ complCumulative d chi2 < p
  | otherwise      = error $ "Statistics.Test.ChiSquare.chi2test: bad p-value: " ++ show p
  where
    n     = G.length vec - ndf - 1
    chi2  = sum $ G.map (\(o,e) -> square (fromIntegral o - e) / e) vec
    d     = chiSquared n
{-# SPECIALIZE
    chi2test :: Double -> Int -> U.Vector (Int,Double) -> TestResult #-}
{-# SPECIALIZE
    chi2test :: Double -> Int -> V.Vector (Int,Double) -> TestResult #-}

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
{-# LANGUAGE OverloadedLists #-}

module Main where

import Numeric.GSL.ODE
import Numeric.LinearAlgebra

-- Differential equation
f :: Double -> [Double] -> [Double]
f t [x, v] = [v, - x + mu * v * (1 - x ^ 2)]

-- Mu scalar, dampening strenth
mu :: Double
mu = 0.1

-- Boundary conditions
ts :: Vector Double
ts = linspace 1000 (0, 50)

-- Use default solver: Embedded Runge-Kutta-Fehlberg (4, 5) method.
vanderpol1 :: [Vector Double]
vanderpol1 = toColumns $ odeSolve f [1, 0] ts

-- Use Runge-Kutta (2,3) solver
vanderpol2 :: [Vector Double]
vanderpol2 = toColumns $ odeSolveV RK2 hi epsAbs epsRel (l2v f) [1, 0] ts
  where
    epsAbs = 1.49012e-08
    epsRel = epsAbs
    hi = (ts ! 1 - ts ! 0) / 100
    l2v f = \t -> fromList . f t . toList

main :: IO ()
main = do
  print vanderpol1
  print vanderpol2

max_stars_repo_name: aavogt/wiwinwlh

max_issues_repo_name: aavogt/wiwinwlh

max_forks_repo_name: aavogt/wiwinwlh


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
module AI.Singularity.Data.Vector where
import GHC.TypeLits
import qualified Numeric.LinearAlgebra.Data as LD
import qualified Numeric.LinearAlgebra.Devel as LDev
import qualified Numeric.LinearAlgebra as LA
import Control.Lens
import Data.Proxy
import Control.Monad.ST
import System.Random
import Control.Monad
import Control.Arrow
import qualified Data.Vector.Storable as VS

import AI.Singularity.Utils
import AI.Singularity.Data.Matrix

data GeneralVector :: (* -> *) -> Nat -> * where
  Vec :: vec Float -> GeneralVector vec n
deriving instance Show (vec Float) => Show (GeneralVector vec n)

vec :: Lens' (GeneralVector vec n) (vec Float)
{-# INLINE vec #-}
vec f (Vec v) = Vec <$> f v

type Vector = GeneralVector LD.Vector

type MutVector s = GeneralVector (LDev.STVector s)

thawV :: Vector n -> ST s (MutVector s n)
thawV = Vec <$< LDev.thawVector . view vec

cmap :: (Float -> Float) -> Vector n -> Vector n
cmap = over vec . LD.cmap

foldV :: KnownNat n => (Float -> Float -> Float) -> Float -> Vector n -> Float
foldV f z = VS.foldr f z . view vec

lenV :: forall n. KnownNat n => Vector n -> Int
lenV _ = natToInt @ n

fromListV :: forall n. KnownNat n => [Float] -> Vector n
fromListV = Vec . LD.fromList . take l . cycle
  where l = natToInt @ n

fromValV :: forall n. KnownNat n => Float -> Vector n
fromValV = fromListV . repeat

toListV :: forall n. KnownNat n => Vector n -> [Float]
toListV (Vec v) = VS.toList v

consV :: forall n. (KnownNat n) => Float -> Vector n -> Vector (n+1)
consV d = Vec . VS.cons d . view vec

tailV :: forall n. (KnownNat n) => Vector (n+1) -> Vector n
tailV = Vec . VS.tail . view vec

takeV :: forall n m. (KnownNat n, KnownNat m, m <= n) => Vector n -> Vector m
takeV = Vec . VS.take m . view vec
  where m = natToInt @ n

dropV :: forall n m. (KnownNat n, KnownNat m, m <= n) => Vector n -> Vector (n - m)
dropV = Vec . VS.drop m . view vec
  where m = natToInt @ m

splitV :: forall k l n. (KnownNat k, KnownNat l, KnownNat n, (k+l) ~ n) => Vector n -> (Vector k, Vector l)
splitV = view vec >>> VS.take k &&& VS.drop k >>> Vec *** Vec
  where k = natToInt @ k

concatVec :: forall n m. (KnownNat n, KnownNat m) => Vector n -> Vector m -> Vector (n + m)
concatVec (Vec v1) (Vec v2) = Vec (v1 VS.++ v2)

(<.>) :: Vector n -> Vector n -> Float
(<.>) (Vec a) (Vec b) = (LA.<.>) a b

(#>) :: Matrix m n -> Vector m -> Vector n
(#>) (Mat m) (Vec v) = Vec (m LA.#> v)

outer :: Vector m -> Vector n -> Matrix n m
outer (Vec a) (Vec b) = Mat (a `LA.outer` b)

appendV :: (KnownNat n, KnownNat m) => Vector n -> Vector m -> Vector (n+m)
appendV (Vec v1) (Vec v2) = Vec $ v1 VS.++ v2

instance forall n. KnownNat n => Num (Vector n) where
  (+) (Vec v1) (Vec v2) = Vec (v1 + v2)
  (*) (Vec v1) (Vec v2) = Vec (v1 * v2)
  negate (Vec v1)       = Vec (LD.cmap negate v1)
  fromInteger i         = fromListV [fromIntegral i]
  abs    (Vec v1)       = Vec (LD.cmap abs v1)
  signum (Vec v1)       = Vec (LD.cmap signum v1)

instance forall n. KnownNat n => Fractional (Vector n) where
  (/) (Vec v1) (Vec v2) = Vec (v1 / v2)
  fromRational          = fromValV . fromRational

instance forall n. KnownNat n => Floating (Vector n) where
  pi    = fromValV pi
  exp   = cmap exp
  log   = cmap log
  sin   = cmap sin
  cos   = cmap cos
  sinh  = cmap sinh
  cosh  = cmap cosh
  tanh  = cmap tanh
  asinh = cmap asinh
  acosh = cmap acosh
  atanh = cmap atanh

answer :: Vector n -> Int
answer (Vec v) = fst . VS.ifoldr (\i x b -> if x > snd b then (i,x) else b) (0,0) $ v

max_stars_repo_name: Antystenes/Memetic-Predictor

max_issues_repo_name: Antystenes/Memetic-Predictor

max_forks_repo_name: Antystenes/Memetic-Predictor


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Data.Random.Distribution.MultivariateNormal
-- Copyright   :  (c) 2016 FP Complete Corporation
-- License     :  MIT (see LICENSE)
-- Maintainer  :  dominic@steinitz.org
--
-- Sample from the multivariate normal distribution with a given
-- vector-valued \(\mu\) and covariance matrix \(\Sigma\). For example,
-- the chart below shows samples from the bivariate normal
-- distribution.
--
-- <<diagrams/src_Data_Random_Distribution_MultivariateNormal_diagM.svg#diagram=diagM&height=600&width=500>>
--
-- Example code to generate the chart:
--
-- > import qualified Graphics.Rendering.Chart as C
-- > import Graphics.Rendering.Chart.Backend.Diagrams
-- >
-- > import Data.Random.Distribution.MultivariateNormal
-- >
-- > import qualified Data.Random as R
-- > import Data.Random.Source.PureMT
-- > import Control.Monad.State
-- > import qualified Numeric.LinearAlgebra.HMatrix as LA
-- >
-- > nSamples :: Int
-- > nSamples = 10000
-- >
-- > sigma1, sigma2, rho :: Double
-- > sigma1 = 3.0
-- > sigma2 = 1.0
-- > rho = 0.5
-- >
-- > singleSample :: R.RVarT (State PureMT) (LA.Vector Double)
-- > singleSample = R.sample $ Normal (LA.fromList [0.0, 0.0])
-- >                (LA.sym $ (2 LA.>< 2) [ sigma1, rho * sigma1 * sigma2
-- >                                      , rho * sigma1 * sigma2, sigma2])
-- >
-- > multiSamples :: [LA.Vector Double]
-- > multiSamples = evalState (replicateM nSamples $ R.sample singleSample) (pureMT 3)
-- > pts = map (f . LA.toList) multiSamples
-- >   where
-- >     f [x, y] = (x, y)
-- >     f _      = error "Only pairs for this chart"
-- >
-- >
-- > chartPoint pointVals n = C.toRenderable layout
-- >   where
-- >
-- >     fitted = C.plot_points_values .~ pointVals
-- >               $ C.plot_points_style  . C.point_color .~ opaque red
-- >               $ C.plot_points_title .~ "Sample"
-- >               $ def
-- >
-- >     layout = C.layout_title .~ "Sampling Bivariate Normal (" ++ (show n) ++ " samples)"
-- >            $ C.layout_y_axis . C.laxis_generate .~ C.scaledAxis def (-3,3)
-- >            $ C.layout_x_axis . C.laxis_generate .~ C.scaledAxis def (-3,3)
-- >
-- >            $ C.layout_plots .~ [C.toPlot fitted]
-- >            $ def
-- >
-- > diagM = do
-- >   denv <- defaultEnv C.vectorAlignmentFns 600 500
-- >   return $ fst $ runBackend denv (C.render (chartPoint pts nSamples) (500, 500))
--
-----------------------------------------------------------------------------

{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE MultiParamTypeClasses #-}

module Data.Random.Distribution.MultivariateNormal
    ( Normal(..)
    ) where

import           Data.Random.Distribution
import qualified Numeric.LinearAlgebra.HMatrix as H
import           Control.Monad
import qualified Data.Random as R
import           Foreign.Storable ( Storable )
import           Data.Maybe ( fromJust )

normalMultivariate :: H.Vector Double -> H.Herm Double -> R.RVarT m (H.Vector Double)
normalMultivariate mu bigSigma = do
  z <- replicateM (H.size mu) (rvarT R.StdNormal)
  return $ mu + bigA H.#> (H.fromList z)
  where
    (vals, bigU) = H.eigSH bigSigma
    lSqrt = H.diag $ H.cmap sqrt vals
    bigA = bigU H.<> lSqrt

data family Normal k :: *

data instance Normal (H.Vector Double) = Normal (H.Vector Double) (H.Herm Double)

instance Distribution Normal (H.Vector Double) where
  rvar (Normal m s) = normalMultivariate m s

normalPdf :: (H.Numeric a, H.Field a, H.Indexable (H.Vector a) a, Num (H.Vector a)) =>
             H.Vector a -> H.Herm a -> H.Vector a -> a
normalPdf mu sigma x = exp $ normalLogPdf mu sigma x

normalLogPdf :: (H.Numeric a, H.Field a, H.Indexable (H.Vector a) a, Num (H.Vector a)) =>
                 H.Vector a -> H.Herm a -> H.Vector a -> a
normalLogPdf mu bigSigma x = - H.sumElements (H.cmap log (diagonals dec))
                              - 0.5 * (fromIntegral (H.size mu)) * log (2 * pi)
                              - 0.5 * s
  where
    dec = fromJust $ H.mbChol bigSigma
    t = fromJust $ H.linearSolve (H.tr dec) (H.asColumn $ x - mu)
    u = H.cmap (\v -> v * v) t
    s = H.sumElements u

diagonals :: (Storable a, H.Element t, H.Indexable (H.Vector t) a) =>
             H.Matrix t -> H.Vector a
diagonals m = H.fromList (map (\i -> m H.! i H.! i) [0..n-1])
  where
    n = max (H.rows m) (H.cols m)

instance PDF Normal (H.Vector Double) where
  pdf (Normal m s) = normalPdf m s
  logPdf (Normal m s) = normalLogPdf m s

max_stars_repo_name: steinitznavican/random-fu-multivariate

max_issues_repo_name: steinitznavican/random-fu-multivariate

max_forks_repo_name: steinitznavican/random-fu-multivariate


text:
module Lib
  (
    trainNetwork
  , classifyTest
  , generateExps
  , runExperiment
  , runExpsIO
  , runExpIO
  ) where

import Numeric.LinearAlgebra

import Net
import CSVMatrixParser

import Data.List
import Text.Printf

type LearningRate = Double
type Epochs = Int
type HiddenUnits = Int
type Sigmoid = Bool

data Experiment = Experiment LearningRate Epochs HiddenUnits Sigmoid deriving (Show, Eq)

latexExp :: Experiment -> String
latexExp (Experiment lr e hu s) = show hu ++ " & " ++ show e ++ " & "
                                  ++ (printf "%.3f" lr :: String) ++ " & " ++ (if s then "Sigmoid" else "Tanh")

latexResults :: [Double] -> String
latexResults results = foldl (\str r -> str ++ " & " ++ (printf "%.2f" (100 * r) :: String) ++ "\\%") "" results

generateExps :: [LearningRate] -> [Epochs] -> [HiddenUnits] -> [Sigmoid] -> [Experiment]
generateExps lrs es hus ss = [Experiment lr e hu s | lr <- lrs, e <- es, hu <- hus, s <- ss]

runExpIO :: String -> [String] -> Experiment -> IO ()
runExpIO trainData testData exp = do
  results1 <- runExperiment trainData testData exp
  results2 <- runExperiment trainData testData exp
  results3 <- runExperiment trainData testData exp
  let avgResults = zipWith3 (\a b c -> (a + b + c) / 3.0) results1 results2 results3
  putStrLn (latexExp exp ++ latexResults avgResults ++ "\\\\")

runExpsIO :: String -> [String] -> [Experiment] -> IO ()
runExpsIO trainData testData = mapM_ (runExpIO trainData testData)

runExperiment :: String -> [String] -> Experiment -> IO [Double]
runExperiment trainData testData (Experiment lr e hu s)  = do
  network <- trainNetwork trainData lr e hu s
  sequence $ map (classifyTest network) testData

trainNetwork :: String -> LearningRate -> Epochs -> HiddenUnits -> Sigmoid -> IO NN
trainNetwork trainingData lr epochs nh isSigmoid = do
  let (f, f') = if isSigmoid then (sigmoid, sigmoid') else (tanh, tanh')
  (i,o) <- getData trainingData
  net   <- buildNetwork (f,f') [5,nh,2]
  return $ repTrain epochs lr i o net

classifyTest :: NN -> String -> IO Double
classifyTest net testData = do
  (i,o) <- getData testData
  let output = forward net i
      testRows = toRows o
      expRows = toRows output
      testClass = map classification testRows
      expClass = map classification expRows
      correct = zipWith xnor testClass expClass
      numCorrect = length $ filter id correct
      percentageCorrect = (fromIntegral numCorrect) / (fromIntegral $ length correct)
  return percentageCorrect

classification :: Vector R -> Bool
classification vec = vec ! 0 < 0.5 && vec ! 1 > 0.5

xnor :: Bool -> Bool -> Bool
xnor True True = True
xnor False False = True
xnor _ _ = False


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
module TensorFlow.Operators
  ( (^+)
  , (^-)
  , (^*)
  , (^/)
  ) where

import Data.Int (Int8, Int16, Int32, Int64)
import Data.Word (Word8, Word16)
import Data.Complex (Complex)
import qualified TensorFlow.Core as TF (Tensor, Build, OneOf)
import qualified TensorFlow.GenOps.Core as TF (add, sub, mul, realDiv)

infixl 6 ^+, ^-
infixl 7 ^*, ^/

(^+) :: TF.OneOf '[Complex Double, Complex Float, Int16, Int32, Int64, Int8, Word16, Word8, Double, Float] t
     => TF.Tensor v'1 t
     -> TF.Tensor v'2 t
     -> TF.Tensor TF.Build t
a ^+ b = a `TF.add` b

(^-) :: TF.OneOf '[Complex Double, Complex Float, Int32, Int64, Word16, Double, Float] t
     => TF.Tensor v'1 t
     -> TF.Tensor v'2 t
     -> TF.Tensor TF.Build t
a ^- b = a `TF.sub` b

(^*) :: TF.OneOf '[Complex Double, Complex Float, Int16, Int32, Int64, Int8, Word16, Word8, Double, Float] t
     => TF.Tensor v'1 t
     -> TF.Tensor v'2 t
     -> TF.Tensor TF.Build t
a ^* b = a `TF.mul` b

(^/) :: TF.OneOf '[Complex Double, Complex Float, Int16, Int32, Int64, Int8, Word16, Word8, Double, Float] t
     => TF.Tensor v'1 t
     -> TF.Tensor v'2 t
     -> TF.Tensor TF.Build t
a ^/ b = a `TF.realDiv` b

max_stars_repo_name: bjoeris/orbits-haskell-tensorflow

max_issues_repo_name: bjoeris/orbits-haskell-tensorflow

max_forks_repo_name: bjoeris/orbits-haskell-tensorflow


text:
module TwoMoon where

  import Numeric.LinearAlgebra
  import System.Environment

  import Neural
  import Test
  import TwoMoon.Load

  main :: IO ()
  main = do
    twomoon <- loadExamples
    learningRate <- (getArgs >>= (\args -> (return.read) $ head args))
    epoch <- (getArgs >>= (\args -> (return.read) $ args!!1))
    let dataSet = splitDataSet twomoon
    runSuite epoch learningRate networkSetting dataSet

  -- two hidden layer with 20, 10 nodes each
  networkSetting = [2, 4]

max_stars_repo_name: ycjungSubhuman/hskmlp

max_issues_repo_name: ycjungSubhuman/hskmlp

max_forks_repo_name: ycjungSubhuman/hskmlp


text:
{-|
Module      : MachineLearning.Model.Regression
Description : TIR expression data structures
Copyright   : (c) Fabricio Olivetti de Franca, 2022
License     : GPL-3
Maintainer  : fabricio.olivetti@gmail.com
Stability   : experimental
Portability : POSIX

Fitting functions for the coefficients.
-}
module MachineLearning.Model.Regression
  ( fitTask
  , predictTask
  , evalPenalty
  , applyMeasures
  , nonlinearFit
  , tirToMatrix
  ) where

import Data.Bifunctor
import           Data.SRTree                   (evalFun, Function(..), OptIntPow(..), evalFun, inverseFunc)
import MachineLearning.TIR       (TIR(..),  Individual(..), Dataset, Constraint, assembleTree, replaceConsts)
import qualified Data.Vector           as V
import qualified Data.Vector.Storable  as VS
import Data.Vector           ((!))
import           Data.List                     (nub)
import           Data.Vector.Storable          (Vector, splitAt)
import           Numeric.LinearAlgebra         ((<\>), Matrix)
import           Numeric.GSL.Fitting           (nlFitting, FittingMethod(..))
import           Prelude               hiding  (splitAt)

import qualified Numeric.LinearAlgebra                     as LA

import MachineLearning.TIR (Individual(..))
import MachineLearning.Utils.Config (Task(..), Penalty(..))
import MachineLearning.Model.Measure (Measure(..))

-- * IT specific stuff

createQ :: Vector Double -> LA.Matrix Double -> LA.Matrix Double
createQ ys = LA.fromColumns . map (*ys) . LA.toColumns
{-# INLINE createQ #-}

avg :: Vector Double -> Vector Double
avg ys = LA.fromList [LA.sumElements ys / fromIntegral (LA.size ys)]
{-# INLINE avg #-}

-- | transform a data matrix using a TIR expression. This function returns
-- a tuple with the transformed data of the numerator and denominator, respectivelly.
-- Each column of the transformed data represents one term of the TIR expression.
tirToMatrix :: Dataset Double -> TIR -> (LA.Matrix Double, LA.Matrix Double)
tirToMatrix xss (TIR _ p q) = bimap (LA.fromColumns . (bias:)) LA.fromColumns (p', q')
  where
    bias      = V.head xss
    xss'      = V.tail xss
    sigma2mtx = map (\(_, g, ps) -> evalFun g $ evalPi ps)
    evalPi    = foldr (\(ix, k) acc -> acc * (xss' ! ix ^^ k)) 1
    p'        = sigma2mtx p
    q'        = sigma2mtx q
{-# INLINE tirToMatrix #-}

-- | Fits a linear model using l2-penalty
ridge :: Matrix Double -> Vector Double -> Matrix Double
ridge a b = oA <\> oB
  where
   mu = 0.01

   a' = LA.tr a
   b' = LA.tr $ LA.asColumn b
   oA = (a' <> LA.tr a') + (mu * LA.ident (LA.rows a'))
   oB = a' <> LA.tr b'
{-# INLINE ridge #-}

-- | Predicts a linear model
predict :: Matrix Double -> Vector Double -> Vector Double
predict xs w | LA.cols xs == LA.size w = xs LA.#> w
             | otherwise = error $ "predict: " ++ show (LA.size xs) ++ show (LA.size w)
{-# INLINE predict #-}

-- | Solve the OLS *zss*w = ys*
solveOLS :: Matrix Double -> Vector Double -> Vector Double
solveOLS zss ys = zss <\> ys
{-# INLINE solveOLS #-}

-- | Applies OLS and returns a Solution
-- if the expression is invalid, it returns Infinity as a fitness
--regress :: Matrix Double -> Vector Double -> [Vector Double]
regress :: TIR -> Dataset Double -> Vector Double -> [Vector Double]
regress tir xss ys = [ws]
  where
    (zssP, zssQ) = tirToMatrix xss tir
    ys'          = evalFun (inverseFunc $ _funY tir) ys
    zssQy        = createQ ys' zssQ
    zss          = if LA.cols zssQ >= 1
                       then zssP LA.||| negate zssQy
                       else zssP
    ws             = if LA.cols zss == 1
                       then avg ys'
                       else solveOLS zss ys'
-- regress zss ys = [solveOLS zss ys]
{-# INLINE regress #-}

-- | Applies conjugate gradient for binary classification
--classify :: Matrix Double -> Vector Double -> [Vector Double]
classify :: Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
classify niter tir xss ys = [ws]
  where
    ws           = nonlinearFit niter zssP zssQ ys sigmoid dsigmoid theta0
    theta0       = LA.konst 0 (LA.cols zssP + LA.cols zssQ)
    (zssP, zssQ) = tirToMatrix xss tir

-- | Applies conjugate gradient for one-vs-all classification
--classifyMult :: Matrix Double -> Vector Double -> [Vector Double]
classifyMult :: Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
classifyMult niter tir xss ys = zipWith minimize yss theta0
  where
    numLabels    = length $ nub $ LA.toList ys
    yss          = map f [0 .. numLabels-1]
    minimize y t = nonlinearFit niter zssP zssQ y sigmoid dsigmoid t
    theta0       = replicate numLabels $ LA.konst 0 (LA.cols zssP + LA.cols zssQ)
    (zssP, zssQ) = tirToMatrix xss tir

    f sample = VS.map (\a -> if round a == sample then 1 else 0) ys

-- | chooses the appropriate fitting function
--fitTask :: Task -> Matrix Double -> Vector Double -> [Vector Double]
fitTask :: Task -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
fitTask Regression             = regress
fitTask (RegressionNL niter)   = regressNL niter
fitTask (Classification niter) = classify niter
fitTask (ClassMult niter)      = classifyMult niter
{-# INLINE fitTask #-}

-- | sigmoid function for classification.
sigmoid :: Floating a => a -> a
sigmoid z = 1 / (1+exp(-z))
{-# INLINE sigmoid #-}

-- | derivative sigmoid function for classification.
dsigmoid :: Floating a => a -> a
dsigmoid z = sigmoid z * (1 - sigmoid z)
{-# INLINE dsigmoid #-}

-- | chooses the appropriate prediction function
predictTask :: Task -> [Vector Double] -> Vector Double
predictTask _ []                   = error "predictTask: empty coefficients matrix"
predictTask Regression yss         = head yss
predictTask (RegressionNL _) yss   = head yss
predictTask (Classification _) yss = sigmoid $ head yss
predictTask (ClassMult _) yss      = LA.vector $ map (fromIntegral . LA.maxIndex) $ LA.toRows $ LA.fromColumns $ map sigmoid yss
{-# INLINE predictTask #-}

-- | evals the penalty function
evalPenalty :: Penalty -> Int -> Double -> Double
evalPenalty NoPenalty _   _   = 0.0
evalPenalty (Len c)   len _   = fromIntegral len * c
evalPenalty (Shape c) _   val = val*c
{-# INLINE evalPenalty #-}

-- | applies a list of performance measures
applyMeasures :: [Measure] -> Vector Double -> Vector Double -> [Double]
applyMeasures measures ys ysHat = map ((`uncurry` (ys, ysHat)) . _fun) measures
{-# INLINE applyMeasures #-}

regressNL :: Int -> TIR -> Dataset Double -> Vector Double -> [Vector Double]
regressNL niter tir xss ys = [ws]
  where
    ws           = nonlinearFit niter zssP zssQ ys f f' theta0
    f            = evalFun $ _funY tir
    f'           = derivative $ _funY tir
    --theta0       = LA.konst 1 (LA.cols zssP + LA.cols zssQ)
    theta0       = head $ regress tir xss ys
    (zssP, zssQ) = tirToMatrix xss tir

-- | Non-linear optimization using Levenberg-Marquardt method.
--nonlinearFit :: Monad m => Vector Double -> Matrix Double -> Matrix Double -> Vector Double -> m (Vector Double)
nonlinearFit :: Int
             -> Matrix Double
             -> Matrix Double
             -> Vector Double
             -> (Vector Double -> Vector Double)
             -> (Vector Double -> Vector Double)
             -> Vector Double
             -> Vector Double
nonlinearFit niter zssP zssQ ys f f' theta0 = fst $ nlFitting LevenbergMarquardtScaled 1e-6 1e-6 niter model' jacob' theta0
  where
    model'       = model f ys zssP zssQ
    jacob'       = jacob f' zssP zssQ

-- | calculates the error given the parameter vector beta
model :: (Vector Double -> Vector Double) -> Vector Double -> Matrix Double -> Matrix Double -> Vector Double -> Vector Double
model f ys zssP zssQ beta
  | LA.cols zssQ == 0 = f ysHat_P - ys
  | otherwise         = f ysHat - ys
  where
    (betaP, betaQ) = splitAt (LA.cols zssP) beta
    ysHat_P        = predict zssP betaP
    ysHat_Q        = if LA.cols zssQ == 0 then 0 else predict zssQ betaQ
    ysHat          = ysHat_P / (1 + ysHat_Q)

-- | calculates the Jacobian given the parameter vector beta. Doesn't support
-- the outer transformation function.
jacob :: (Vector Double -> Vector Double) -> Matrix Double -> Matrix Double -> Vector Double -> Matrix Double
jacob f zssP zssQ beta | LA.cols zssQ == 0 = zssP
                       | otherwise         = LA.fromColumns $ pjac <> qjac
  where
    (betaP, betaQ) = splitAt (LA.cols zssP) beta
    ysHat_P        = predict zssP betaP
    ysHat_Q        = predict zssQ betaQ
    ysHat          = f $ ysHat_P / (1 + ysHat_Q)
    pjac           = [ysHat * c / ysHat_Q | c <- LA.toColumns zssP]
    qjac           = [-ysHat * c * (ysHat_P / (1 + ysHat_Q)^2) | c <- LA.toColumns zssQ]

derivative :: (Eq val, Floating val) => Function -> val -> val
derivative Id      = const 1
derivative Abs     = \x -> x / abs x
derivative Sin     = cos
derivative Cos     = negate.sin
derivative Tan     = recip . (**2.0) . cos
derivative Sinh    = cosh
derivative Cosh    = sinh
derivative Tanh    = (1-) . (**2.0) . tanh
derivative ASin    = recip . sqrt . (1-) . (^2)
derivative ACos    = negate . recip . sqrt . (1-) . (^2)
derivative ATan    = recip . (1+) . (^2)
derivative ASinh   = recip . sqrt . (1+) . (^2)
derivative ACosh   = \x -> 1 / (sqrt (x-1) * sqrt (x+1))
derivative ATanh   = recip . (1-) . (^2)
derivative Sqrt    = recip . (2*) . sqrt
derivative Square  = (2*)
derivative Exp     = exp
derivative Log     = recip
{-# INLINE derivative #-}
-- (w1 * p1 + w2 * p2) / (1 + w3 * p3 + w4 * p4)
-- d/dw1 = p1 / (1 + w3 * p3 + w4 * p4)
-- d/dw2 = p2 / (1 + w3 * p3 + w4 * p4)
-- d/dw3 = -p3 * (w1 * p2 + w2 * p2) / (1 + w3 * p3 + w4 * p4)^2
{-
toEv :: SRTree Int Double -> (V.Vector Double -> Double)
toEv (Var !ix) = (`V.unsafeIndex` ix)
toEv (Const !val) = const val
toEv (Add !l !r) = jn (+) (toEv l) (toEv r)
toEv (Mul !l !r) = jn (*) (toEv l) (toEv r)
toEv (Fun Exp !t) = exp . toEv t
{-# INLINE toEv #-}

jn :: (Double -> Double -> Double) -> (V.Vector Double -> Double) -> (V.Vector Double -> Double) -> (V.Vector Double -> Double)
jn op f g = \x -> op (f x) (g x)
{-# INLINE jn #-}
-}

max_stars_repo_name: folivetti/tir

max_issues_repo_name: folivetti/tir

max_forks_repo_name: folivetti/tir


text:
{-# LANGUAGE OverloadedStrings #-}

module Haskem.Lib
    ( doGaussianParsing,
      doFreqParsing,
      doGenericParsing,
      xyzParser,
      gaussianFreqParser,
      gaussianParser,
      elemMass,
      elemSym2Num,
      calculateCoM,
      calcDisplacementVector,
      writeMoleculeToXYZ,
      projectDonN,
      calcProjectionVec,
      calcLength3D,
      setUpBonds,
      setUpBonds',
      setUpBonds'',
      setUpBonds''',
      checkBondCompleteness
    ) where

import              Haskem.Parsers
import              Haskem.Types
import              Control.Applicative
import qualified    Formatting                  as F
import              Data.Attoparsec.Text.Lazy   hiding (take)
import qualified    Data.Text.Lazy              as T
import qualified    Data.Text.IO                as TIO
import              Data.Char
import              Lens.Micro.Platform
import              Debug.Trace
import              Data.List                   hiding (takeWhile)
import              Data.List.Extra
-- import              Prelude                     hiding (takeWhile)
-- import              Data.Tuple.HT               hiding (double, triple)
import qualified    Data.Map.Strict             as Map
import qualified    Data.Vector.Storable        as VS
import qualified    Numeric.LinearAlgebra.Data  as LA
import qualified    Data.IntMap                 as IM
import qualified    Data.IntSet                 as IS

xyzFormat :: F.Format r (String -> Double -> Double -> Double -> r)
xyzFormat = F.string F.% "\t" F.% F.fixed 4 F.% "\t" F.% F.fixed 4 F.% "\t" F.% F.fixed 4

doGaussianParsing :: FilePath -> IO GaussianInfo
doGaussianParsing filePath = do
    fileContent <- TIO.readFile filePath
    let gaussParse = parseOnly gaussianParser fileContent

    case gaussParse of
        Left    err     -> error err
        Right   gInfo   -> return gInfo

-- Parses a Gaussian freq calculation output to a GaussianInfo instance and a
-- list of Normalmode instances. Returns the output as a tuple in IO.
doFreqParsing :: FilePath -> IO (GaussianInfo, [NormalMode])
doFreqParsing filePath = do
    fileContent <- TIO.readFile filePath

    let gaussParse = parseOnly gaussianParser       fileContent
    let gaussFreq  = parseOnly gaussianFreqParser   fileContent

    gFreq <- case gaussFreq of
        Left    err     -> error err
        Right   gFreq   -> return gFreq

    gInfo <- case gaussParse of
        Left    err     -> error err
        Right   gInfo   -> return gInfo

    return (gInfo, gFreq)

-- Parses a file using a specified parser and returns the result in IO
doGenericParsing :: FilePath -> Parser a -> IO a
doGenericParsing filePath parser = do
    fileContent <- TIO.readFile filePath
    let gParse = parseOnly parser fileContent

    case gParse of
        Left    err     -> error err
        Right   gInfo   -> return gInfo

getListOfCoords :: Molecule -> [CartesianCoord]
getListOfCoords molA' = molA' ^.. atoms . each . coordinate

-- Calculates the center of mass of a given molecule of type Molecule (from Haskem.Types)molA
calculateCoM :: Molecule -> CartesianCoord
calculateCoM molecule' = CartesianCoord xCoM yCoM zCoM
    where
        atomsList   = getListOfCoords molecule'
        massList    = map elemMass $ molecule' ^.. atoms . each . atomNumber
        fullMass    = sum massList
        xCoM        = (/ fullMass) $ sum $ zipWith (*) massList $ map (^. xCoord) atomsList
        yCoM        = (/ fullMass) $ sum $ zipWith (*) massList $ map (^. yCoord) atomsList
        zCoM        = (/ fullMass) $ sum $ zipWith (*) massList $ map (^. zCoord) atomsList

-- Calculate a displacement vector from two molecular geometries.
-- Outputs a maybe molecule (in case the geometries do not match
-- (checked via number of atoms and the list of element numbers))
calcDisplacementVector :: Molecule -> Molecule -> [CartesianCoord]
calcDisplacementVector molA molB =
    do
        let nrOfAtomsA  = molA ^. nrOfAtoms
            nrOfAtomsB  = molA ^. nrOfAtoms
            elemA   = molA ^.. atoms . each . atomNumber
            elemB   = molB ^.. atoms . each . atomNumber
        if (nrOfAtomsA == nrOfAtomsB) && (elemA == elemB)
            then do
                let atomsA  = getListOfCoords molA
                    atomsB  = getListOfCoords molB
                    xD      = zipWith (-) (map (^. xCoord) atomsA) $ map (^. xCoord) atomsB
                    yD      = zipWith (-) (map (^. yCoord) atomsA) $ map (^. yCoord) atomsB
                    zD      = zipWith (-) (map (^. zCoord) atomsA) $ map (^. zCoord) atomsB
                zipWith3 CartesianCoord xD yD zD
            else
                []

-- Writes a "Molecule" back to disk as a .xyz file
writeMoleculeToXYZ :: Molecule -> FilePath -> IO ()
writeMoleculeToXYZ molA fp = writeFile fp (show nrOfAtomsA
                                ++ "\nThanks for using Haskem!\n"
                                ++ T.unpack (T.unlines lineList))
    where
        coordinates = getListOfCoords molA
        nrOfAtomsA  = molA ^. nrOfAtoms
        elems       = map elemNum2Sym $ molA ^.. atoms . each . atomNumber
        lineList    = zipWith4 (F.format xyzFormat)
                        elems
                        (map (^. xCoord) coordinates)
                        (map (^. yCoord) coordinates)
                        (map (^. zCoord) coordinates)

projectDonN :: [CartesianCoord] -> NormalMode -> Double
projectDonN displCoordsA modesA = projectCoordNM
    where
        displCoordsNM   = modesA ^.. displacement . each .coordinate  :: [CartesianCoord]
        projectCoordNM  = sum [sum $ zipWith (*)
                            (map ( ^. cc) displCoordsA)
                            (map ( ^. cc) displCoordsNM) | cc <- [xCoord, yCoord, zCoord]]

calcProjectionVec :: [NormalMode] -> [CartesianCoord] -> [Double]
calcProjectionVec nModes displCoords = projectionVec
        where
            projectionVec = map (projectDonN displCoords) nModes

{-
Build redundant internal coordinates from the molecular cartesian structure
-}

calcLength3D :: CartesianCoord -> CartesianCoord -> Double
calcLength3D coord1 coord2 = length3D
        where
            length3D = sqrt . sum $ (** 2) <$> dVector
            dVector = [(coord1 ^. axis) - (coord2 ^. axis)
                            | axis <- [xCoord, yCoord, zCoord]]

setUpBonds :: Molecule -> [(String, String, Double)]
setUpBonds molA = bondsList
        where
            atomListN   = molA ^.. atoms . each . atomNumber
            atomListS   = map elemNum2Sym   atomListN
            vdWList     = map covRadius     atomListN
            enumCoords  =  zip4 [1..] atomListS vdWList $ getListOfCoords molA
            bondsList   = [(e1 ++ show i1, e2 ++ show i2, calcLength3D c1 c2)
                            | (i1, e1, r1, c1) <- enumCoords,
                              (i2, e2, r2, c2) <- enumCoords,
                              let l12 = calcLength3D c1 c2,
                              i1 < i2 && l12 < 1.3 * (r1 + r2)
                            ]

setUpBonds' :: Molecule -> [LinearBond]
setUpBonds' molA = bondsList
        where
            atomListN'   = molA ^.. atoms . each . atomNumber
            atomListS   = map elemNum2Sym   atomListN'
            vdWList     = map covRadius     atomListN'
            enumCoords  =  zip4 [1..] atomListS vdWList $ getListOfCoords molA
            bondsList   = [LinearBond (AtomSpec e1 i1) (AtomSpec e2 i2) l12
                            | (i1, e1, r1, c1) <- enumCoords,
                              (i2, e2, r2, c2) <- enumCoords,
                              i1 < i2,
                              let l12 = calcLength3D c1 c2,
                              l12 <= 1.3 * (r1+r2)]

setUpBonds'' :: Molecule -> LA.Matrix Double
setUpBonds'' molA = LA.reshape nrOfAtoms' bondsMatrix
        where
            nrOfAtoms'  = molA ^. nrOfAtoms
            atomListN'  = molA ^.. atoms . each . atomNumber
            atomListS   = map elemNum2Sym   atomListN'
            vdWList     = map covRadius     atomListN'
            enumCoords  = zip4 [1..] atomListS vdWList $ getListOfCoords molA
            bondsMatrix = VS.fromList [bl | (i1, e1, r1, c1) <- enumCoords,
                                            (i2, e2, r2, c2) <- enumCoords,
                                            let l12 = calcLength3D c1 c2,
                                            let bl  = if l12 < 1.3 * (r1+r2) then 1 else 0]

setUpBonds''' :: Molecule -> [IS.IntSet]
setUpBonds''' molA = bondsMatrix
        where
            nrOfAtoms'  = molA ^. nrOfAtoms
            atomListN'  = molA ^.. atoms . each . atomNumber
            atomListS   = map elemNum2Sym   atomListN'
            vdWList     = map covRadius     atomListN'
            enumCoords  = zip4 [1..] atomListS vdWList $ getListOfCoords molA
            bondsMatrix = IS.fromList <$> chunksOf nrOfAtoms' [bl |
                                                (i1, e1, r1, c1) <- enumCoords,
                                                (i2, e2, r2, c2) <- enumCoords,
                                                let l12 = calcLength3D c1 c2,
                                                let bl  = if l12 < 1.3 * (r1+r2) then 1 else 0]
-- TODO!

removeDuplicates :: (Ord a) => [a] -> [a]
removeDuplicates = map head . group . sort

checkBondCompleteness :: Molecule -> [LinearBond] -> [Atom]
checkBondCompleteness mol' listOfBonds' = missing
        where
            nrOfAtoms'      = mol' ^. nrOfAtoms
            atomNums'       = removeDuplicates $ map (^. bAtomA . geomNum) listOfBonds'
                                ++ map (^. bAtomB . geomNum) listOfBonds'
            missingIndxs    = findIndices (`notElem` atomNums')  [1..nrOfAtoms']
            missing         = case missingIndxs of
                               []   -> []
                               _    -> concat [mol' ^.. atoms . ix ii | ii <- missingIndxs]

defineFragments     :: [IS.IntSet] -> [IS.IntSet]
defineFragments []      = []
defineFragments [x]     = [x]

defineFragments (a:as)  =
    if (all null overlaps)
        then (a:as)
        else defineFragments reducedSet
    where
        setsWithOverlap ::  IS.IntSet -> [IS.IntSet] -> [IS.IntSet]
        setsWithOverlap bond bondsList = [ii | ii <- bondsList,
                                                not . IS.null $ IS.intersection bond ii,
                                                bond /= ii]

        overlaps =
            [setsWithOverlap ((a:as) !! ii) (deleteNth ii (a:as))
                | ii <- [0..length(a:as) - 1]
            ]
            :: [[IS.IntSet]]

        reducedSet =
            nub
            [ IS.unions (((a:as) !! i) : (overlaps !! i))
                | i <- [0 .. length (a:as) - 1]
            ]

deleteNth :: Int -> [a] -> [a]
deleteNth n l = (take n l) ++ (drop (n + 1) l)

-- defineFragments'    :: VS.Vector (IM.IntMap Int) -> VS.Vector (IM.IntMap Int)
-- defineFragments'

max_stars_repo_name: FabianGD/Haskem

max_issues_repo_name: FabianGD/Haskem

max_forks_repo_name: FabianGD/Haskem