{-# OPTIONS_GHC -fno-warn-type-defaults #-}
{-# LANGUAGE DataKinds, FlexibleContexts, FlexibleInstances, GADTs #-}
{-# LANGUAGE MultiParamTypeClasses, NoMonomorphismRestriction      #-}
{-# LANGUAGE OverloadedStrings, PolyKinds, TypeFamilies            #-}
{-# LANGUAGE UndecidableInstances                                  #-}
{-# OPTIONS_GHC -fno-warn-type-defaults -fno-warn-orphans -Wall #-}
module Main where
import           Algebra.Bridge.Singular
import           Algebra.Prelude.Core
import           Control.Monad                   (void)
import           Control.Parallel.Strategies
import qualified Data.Text                       as T
import qualified Data.Vector.Unboxed             as V
import           Numeric.Field.Fraction          (Fraction)
import           Prelude                         (read)
import           Statistics.Resampling
import           Statistics.Resampling.Bootstrap
import           Statistics.Types
import qualified System.Random.MWC               as Rand


i2 :: [OrderedPolynomial (Fraction Integer) Grevlex 5]
i2 =  [35 * y^4 - 30*x*y^2 - 210*y^2*z + 3*x^2 + 30*x*z - 105*z^2 +140*y*t - 21*u
      ,5*x*y^3 - 140*y^3*z - 3*x^2*y + 45*x*y*z - 420*y*z^2 + 210*y^2*t -25*x*t + 70*z*t + 126*y*u
      ]
     where [t,u,x,y,z] = vars

i3 :: [OrderedPolynomial (Fraction Integer) Grevlex 4]
i3 = [ x^31 - x^6 - x- y, x^8 - z, x^10 -t]
  where
    [t,x,y,z] = vars

i4 :: [OrderedPolynomial (Fraction Integer) Grevlex 4]
i4 = [ w+x+y+z, w*x+x*y+y*z+z*w, w*x*y + x*y*z + y*z*w + z*w*x, w*x*y*z - 1]
  where
    [x,y,z,w] = vars

benchIdeal :: IsSingularPolynomial poly => Text -> Ideal poly -> IO Double
benchIdeal fun i = fmap ((/1000) . read . T.unpack) $ singular $
  toProg fun i

toProg :: (SingularOrder (Arity poly) (MOrder poly), SingularCoeff (Coefficient poly), IsOrderedPolynomial poly) => Text -> Ideal poly -> Text
toProg fun i =
  let expr = funE fun [idealE' i]
  in prettySingular $ do
    directC "system(\"--ticks-per-sec\", 1000000)"
    void $ ringC "R" expr
    declOnlyC IdealT "G"
    directC "timer=1"
    directC "int t=rtimer"
    letC "G" expr
    directC "print(rtimer-t)"
    directC "exit"

analyse :: (IsSingularPolynomial poly) => String -> Text -> Ideal poly -> IO ()
analyse lab fun ideal = do
  gen <- Rand.create
  i2Gr <- V.replicateM 200 $ benchIdeal fun ideal
  res  <- resample gen [Mean, StdDev] 1000 i2Gr
  let [Estimate mn (ConfInt lbmn ubmn _) ,Estimate dv _]
        = bootstrapBCA cl95 i2Gr res
  putStrLn lab
  putStrLn $ unlines $ map ('\t':)
    ["Mean:\t" ++ show mn ++ "(ms)"
    ,"MeanLB:\t" ++ show lbmn
    ,"MeanUB:\t" ++ show ubmn
    ,"StdDev:\t" ++ show dv
    ]

runTestCases :: (IsMonomialOrder n o, SingularCoeff k, KnownNat n)
             => String -> Ideal (OrderedPolynomial k o n) -> IO ()
runTestCases lab i = do
  analyse (lab ++ " (Lex, Sing(gr))") "groebner" $ fmap (changeOrder Lex) i
  analyse (lab ++ " (Lex, Sing(sba))") "sba" $ fmap (changeOrder Lex) i
  analyse (lab ++ " (Grevlex, Sing(gr))") "groebner" $ fmap (changeOrder Grevlex) i
  analyse (lab ++ " (Grevlex, Sing(sba))") "sba" $ fmap (changeOrder Grevlex) i

main :: IO ()
main = do
  ideal2 <- return $! (toIdeal i2 `using` rdeepseq)
  ideal3 <- return $! (toIdeal i3 `using` rdeepseq)
  ideal4 <- return $! (toIdeal i4 `using` rdeepseq)

  runTestCases "I1" ideal2
  runTestCases "Cyclic4" ideal4
  analyse "I3 (Grevlex, Sing(gr))" "groebner" ideal3
  analyse "I3 (Grevlex, Sing(sba))" "sba" ideal3
  -- analyse "groebner" ideal4
  -- analyse "sba" ideal4
  -- analyse "groebner" ideal3
  -- analyse "sba" ideal3
  -- defaultMainWith defaultConfig $
  --      mkTestCases (Left 1) ideal2
  --   ++ mkTestCases (Right "Cyclic-4") ideal4
  --   ++ [mkTC "grevlex03" ideal3]
