Main.hs

{-# LANGUAGE ScopedTypeVariables #-}

module Main where

-- base
import Control.Arrow
  ( second )
import Control.Monad
  ( when )
import Data.Foldable
  ( for_ )
import qualified Data.List.NonEmpty as NE
  ( NonEmpty(..)
  , fromList, head, length, sort
  )
import Data.Semigroup
  ( Arg(..), Sum(..), Product(..) )
import Data.Traversable
  ( for )
import GHC.Clock
  ( getMonotonicTime )
import Numeric
  ( showFFloat )

-- code-page
import System.IO.CodePage
  ( withCP65001 )

-- containers
import qualified Data.IntMap.Strict as IntMap
  ( fromList, toList )
import Data.Tree
  ( foldTree )

-- deepseq
import Control.DeepSeq
  ( rnf )

-- brush-strokes
import Math.Bezier.Stroke
import Math.Interval
import Math.Linear
import Math.Root.Isolation
import Math.Root.Isolation.Core
import Math.Root.Isolation.Newton
import Math.Root.Isolation.Newton.GaussSeidel

-- brush-strokes bench:cusps
import Bench.Cases
import Bench.Types

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

pattern ROUND_UP :: Word
pattern ROUND_UP = 0x4000

foreign import ccall "set_sse_rounding_mode"
  setSSERoundingMode :: Word -> IO ()

main :: IO ()
main = withCP65001 $ do
  setSSERoundingMode ROUND_UP
  putStrLn "Running cusp-finding benchmarks."
  for_ benchGroups $ \ ( groupName, benchGroupCases ) -> do
    putStrLn $ unlines
      [ replicate 40 '='
      , "Benchmark group '" ++ groupName ++ "':" ]
    testsWithTime <- for benchGroupCases $ \ tst -> do
      dt <- benchTestCase groupName tst
      return $ Arg dt tst
    let Arg bestTime bestTest = NE.head $ NE.sort testsWithTime
    when ( NE.length benchGroupCases >= 1 ) $
      putStrLn $ unlines $
        [ "Best time in '" ++ groupName ++ "' group: " ++ show bestTime ++ "s" ]
        ++ [ "  (" ++ descr ++ ")"
           | let descr = testDescription bestTest
           , not ( null descr )
           ]

-- TODO: run the benchmark multiple times to reduce noise.
benchTestCase :: String -> TestCase -> IO Double
benchTestCase testName ( TestCase { testDescription, testBrushStroke, testCuspOptions, testStartBoxes } ) = do
  putStr $ unlines
    [ "  " ++ replicate 38 '-'
    , "  -- Test case: " ++ testName ++ ( if null testDescription then "" else " (" ++ testDescription ++ ")" )
    , "  --" ]
  before <- getMonotonicTime
  let ( _, cornerAndCuspFnI ) = brushStrokeFunctions testBrushStroke
      ( trees, dunno, sols ) =
        foldMap
          ( \ ( i, ( trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = mbCusps } ) ) ->
              ( map ( i, ) trees, map ( ( i , ) . fst ) mbCusps, map ( i, ) defCusps ) ) $
          IntMap.toList $
            findCuspsIn testCuspOptions ( snd . cornerAndCuspFnI ) $
              IntMap.fromList
                [ ( i, boxes )
                | ( i, boxes ) <- testStartBoxes
                ]
  rnf dunno `seq` rnf sols `seq` return ()
  after <- getMonotonicTime
  let dt = after - before
  putStrLn $ unlines $
    showResultsList "sols" sols
    ++
    showResultsList "dunno" dunno
    ++
    [ "  --   Time elapsed: " ++ show dt ++ "s"
    , "  --   Tree size: " ++ show (getSum $ foldMap (foldMap sizeRootIsolationTree) trees)
    , "  --   Newton stats: "
    ] ++ map ( "  --     " ++ ) ( showNewtonStats (foldMap (foldMap newtonStats) trees) )
  return dt

-- | Format a list of results with empty set symbol for empty lists
showResultsList :: Show a => String -> [a] -> [String]
showResultsList label items =
  [ "  --   " ++ label ++ ":" ++ (if null items then " ∅" else "") ]
  ++
  [ "  --     • " ++ show item | item <- items ]

sizeRootIsolationTree :: ( Box 2, RootIsolationTree ( Box 2 ) ) -> Sum Int
sizeRootIsolationTree ( box, tree ) = foldMap ( const $ Sum 1 ) ( showRootIsolationTree box tree )


data NewtonStats
  = NewtonStats
  { newtonImprovement :: !( Sum Double )
  , newtonElimination :: !( Sum Int, BigBoxes )
  , newtonTime        :: !( Sum Double )
  , newtonTotal       :: !( Sum Int )
  }
instance Semigroup NewtonStats where
  NewtonStats imp1 el1 t1 tot1 <> NewtonStats imp2 el2 t2 tot2 =
    NewtonStats ( imp1 <> imp2 ) ( el1 <> el2 ) ( t1 <> t2 ) ( tot1 <> tot2 )
instance Monoid NewtonStats where
  mempty = NewtonStats mempty mempty mempty mempty

showNewtonStats :: NewtonStats -> [ String ]
showNewtonStats ( NewtonStats ( Sum improv ) ( Sum elims, big ) ( Sum time ) ( Sum tot ) )
  | tot == 0
  = [ ]
  | otherwise
  = [ "average improvement:  " ++ showPercent ( improv / fromIntegral tot )
    , "average eliminations: " ++ showPercent ( fromIntegral elims / fromIntegral tot )
    , "time elapsed: " ++ show ( time / fromIntegral tot ) ++ "s"
    ]

showPercent :: Double -> String
showPercent x = fixed 3 2 ( 100 * x ) <> "%"

fixed :: Int -> Int -> Double -> String
fixed digitsBefore digitsAfter x =
  case second ( drop 1 ) . break ( == '.' ) $ showFFloat ( Just digitsAfter ) x "" of
    ( as, bs ) ->
      let
        l, r :: Int
        l = length as
        r = length bs
      in
        replicate ( digitsBefore - l ) ' ' <> as <> "." <> bs <> replicate ( digitsAfter - r ) '0'

data BigBoxes =
  BigBoxes
    { biggestArea :: Maybe ( Box 2 )
    , biggestX :: Maybe ( Box 2 )
    , biggestY :: Maybe ( Box 2 )
    }
  deriving stock ( Show, Eq )
instance Semigroup BigBoxes where
  BigBoxes area1 x1 y1 <> BigBoxes area2 x2 y2 =
    BigBoxes
      ( pickBig boxArea area1 area2 )
      ( pickBig widthX x1 x2 )
      ( pickBig widthY y1 y2 )
    where
      pickBig _    Nothing      Nothing   = Nothing
      pickBig _    (Just x1)    Nothing   = Just x1
      pickBig _    Nothing      (Just x2) = Just x2
      pickBig f j1@(Just x1) j2@(Just x2)
        | f x1 >= f x2
        = j1
        | otherwise
        = j2
      widthX ( 𝕀ℝ2 x  _y ) = width x
      widthY ( 𝕀ℝ2 _x y  ) = width y
instance Monoid BigBoxes where
  mempty = BigBoxes Nothing Nothing Nothing


newtonStats :: ( Box 2, RootIsolationTree ( Box 2 ) ) -> NewtonStats
newtonStats ( _box, RootIsolationLeaf {} ) = mempty
newtonStats (  box, RootIsolationStep step boxes )
  | IsolationStep @Newton ( _, TimeInterval dt ) <- step
  = thisImprovement dt <> foldMap newtonStats boxes
  | otherwise
  = foldMap newtonStats boxes
  where
    thisImprovement dt =
      NewtonStats
        { newtonImprovement = Sum $ ( old - new ) / old
        , newtonElimination =
            if new == 0
            then ( Sum 1, BigBoxes ( Just box ) ( Just box ) ( Just box ) )
            else mempty
        , newtonTime = Sum dt
        , newtonTotal = Sum 1
        }
      where
        old = boxArea box
        new = sum ( map ( boxArea . fst ) boxes )

benchGroups :: [ ( String, NE.NonEmpty TestCase ) ]
benchGroups =
  [ ( "ellipse"
    , NE.fromList
      [ ( ellipseTestCase ( mkOpts gs )
          ( 0, 1 ) pi
          ( defaultStartBoxes [ 0 .. 3 ] )
        ) { testDescription = show gs }
      | gs <- [ GS_Partial, GS_Complete ]
      ]
    )
  , ( "trickyCusp2", NE.fromList [ trickyCusp2TestCase ] )
  , ( "rotation", NE.fromList [ ( rotationTestCase ( mkOpts gs ) ) { testDescription = show gs }
                              | gs <- [ GS_Partial, GS_Complete ] ] )
  , ( "missingCusp2", NE.fromList [ missingCusp2TestCase ] )
  ]
  where
    minWidth = 1e-5
    ε_bis = 5e-3
    mkOpts gs_opts =
      RootIsolationOptions
        { rootIsolationAlgorithms = \ hist ->
            defaultRootIsolationAlgorithms minWidth ε_bis
                          ( newtOpts gs_opts hist ) hist
        }
    newtOpts gs_opts hist =
      NewtonLP
      --NewtonGaussSeidel $
      --  ( defaultGaussSeidelOptions @2 @3 hist )
      --    { gsUpdate = gs_opts }

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

{- -- Old testing code.

getR1 (ℝ1 u) = u

eval
  :: ( I i 1 -> Seq ( I i 1 -> StrokeDatum k i ) )
  -> ( I i 1, Int, I i 1 )
  -> StrokeDatum k i
eval f ( t, i, s ) = ( f t `Seq.index` i ) s

mkVal :: Double -> Int -> Double -> ( ℝ 1, Int, ℝ 1 )
mkVal t i s = ( ℝ1 t, i, ℝ1 s )

mkI :: ( Double, Double ) -> 𝕀
mkI ( lo, hi ) = 𝕀 lo hi

mkBox :: ( Double, Double ) -> ( Double, Double ) -> Box 2
mkBox ( t_min, t_max ) ( s_min, s_max ) =
  ( 𝕀 ( ℝ2 t_min s_min ) ( ℝ2 t_max s_max ) )

potentialCusp :: StrokeDatum 3 𝕀 -> Bool
potentialCusp
  ( StrokeDatum
    { ee = D22 { _D22_v = 𝕀 ( ℝ1 ee_min ) ( ℝ1 ee_max ) }
    , 𝛿E𝛿sdcdt = D12 { _D12_v = T ( 𝕀 ( ℝ2 vx_min vy_min ) ( ℝ2 vx_max vy_max ) )}
    }
  ) =  ee_min <= 0 && ee_max >= 0
    && vx_min <= 0 && vx_max >= 0
    && vy_min <= 0 && vy_max >= 0

dEdsdcdt :: StrokeDatum k i -> D ( k - 2 ) ( I i 2 ) ( T ( I i 2 ) )
dEdsdcdt ( StrokeDatum { 𝛿E𝛿sdcdt = v } ) = v

width :: 𝕀ℝ 1 -> Double
width (𝕀 (ℝ1 lo) (ℝ1 hi)) = hi - lo

negV2 :: 𝕀ℝ 2 -> 𝕀ℝ 2
negV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
  let !( 𝕀 x'_lo x'_hi ) = negate $ 𝕀 x_lo x_hi
      !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi
  in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi )

midV2 :: 𝕀ℝ 2 -> ℝ 2
midV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) =
  ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) )


logLines :: [ String ]
logLines =
  [ "E, dE/ds * dc/dt"
  , "{" ++
    (intercalate ","
     [ "{" ++ showD t ++ "," ++ showD s ++ ",{" ++ intercalate "," vals ++ "}}"
     | t <- [ 0.5484, 0.5484 + 0.00001 .. 0.5488 ]
     , s <- [ 0.5479, 0.5479 + 0.00001 .. 0.5483 ]
     , let StrokeDatum
             { ee = D22 ee _ _ _ _ _
             , 𝛿E𝛿sdcdt = D12 (T f) (T (T f_t)) (T (T f_s))
             } = (curvesI (singleton (ℝ1 t)) `Seq.index` i) (singleton (ℝ1 s))
           ℝ2 vx vy = midPoint2 f
           --ℝ2 vx_t vy_t = midPoint2 f_t
           --ℝ2 vx_s vy_s = midPoint2 f_s
           vals = [ showD ( midPoint ee )
                  , "{" ++ showD vx ++ "," ++ showD vy ++ "}"
                --  , "{" ++ showD vx_t ++ "," ++ showD vy_t ++ "}"
                --  , "{" ++ showD vx_s ++ "," ++ showD vy_s ++ "}"
                  ]
     ]
    ) ++ "}"
  ]
  where
    i = 2
    ( curves, curvesI ) = brushStrokeFunctions $ ellipseBrushStroke ( 0, 1 ) pi
    midPoint (𝕀 (ℝ1 lo) (ℝ1 hi)) = 0.5 * ( lo + hi )
    midPoint2 (𝕀 (ℝ2 lo_x lo_y) (ℝ2 hi_x hi_y))
      = ℝ2 ( 0.5 * ( lo_x + hi_x ) ) ( 0.5 * ( lo_y + hi_y ) )

_x ( ℝ2 x _ ) = x
_y ( ℝ2 _ y ) = y

showD :: Double -> String
showD float = showFFloat (Just 6) float ""

-}