content: {-# LANGUAGE ViewPatterns #-}
  {-# OPTIONS_GHC -fno-warn-type-defaults #-}
  {-# LANGUAGE CPP #-}
  
  #if __GLASGOW_HASKELL__ == 810
  -- Work around to fix GHC Issue #15304, issue popped up again in GHC 8.10, it should be fixed in GHC 8.12
  -- This code is meant to reproduce MR 2608 for GHC 8.10
  {-# OPTIONS_GHC -funfolding-keeness-factor=1 -funfolding-use-threshold=80 #-}
  #endif
  
  -------------------------------------------------------------------
  -- |
  -- Copyright   :  (c) 2017-2020 Nathan Waivio
  -- License     :  BSD3
  -- Maintainer  :  Nathan Waivio <nathan.waivio@gmail.com>
  -- 
  -- A program to test Algebra.Geometric.Cl3
  -- The code runs tests on some standard test input and then
  -- runs quckcheck for some trig identities.
  -- 
  -------------------------------------------------------------------
  
  module Main (main) where
  
  import Algebra.Geometric.Cl3
  import Control.Monad (replicateM)
  import Criterion.Main (defaultMain, bench, nfIO, env, Benchmark)
  import System.Random (randomRIO)
  
  ------------------------------------------------------------------
  -- |
  -- This program verifies the approximate equality of various trig
  -- identities to the with the following limitations:
  -- 
  -- * The magnitude of the cliffor is limited in some cases.
  --
  -- * The imaginary part of the eigenvalues are unwrapped, due to the cyclical nature of some of the results, in a few cases.
  --
  -- * The poles of the functions are excluded.
  --
  -- * The poles of the derivatives of the functions are excluded when the cliffor is has a nilpotent component.
  --
  -- * Approximate equivalence is tested due to limitations with respect to floating point math.
  --
  -- 
  -- The following properties are verified in this module:
  --
  -- * log.exp Identity
  --
  -- * exp.log Identity
  --
  -- * abs*signum law
  --
  -- * The definition of recip
  --
  -- * recip.recip Identity
  --
  -- * sin.asin Identity
  --
  -- * asin.sin Identity
  --
  -- * cos.acos Identity
  --
  -- * acos.cos Identity
  --
  -- * sinh.asinh Identity
  --
  -- * asinh.sinh Identity
  --
  -- * cosh.acosh Identity
  --
  -- * acosh.cosh Identity
  --
  -- * Double Sin Identity
  --
  -- * Double Cos Identity
  --
  -- * Double Tan Identity
  --
  -- * Double Sinh Identity
  --
  -- * Double Cosh Identity
  --
  -- * Double Tanh Identity
  --
  -- * Positive Sin Shift Identity
  --
  -- * Negative Sin Shift Identity
  --
  -- * sin^2+cos^2 Identity
  --
  -- * cosh^2-sinh^2 Identity
  --
  -- * Symmetry of Cosh
  --
  -- * Symmetry of Sinh
  --
  -- * Double I Sin
  --
  -- * Composition Algebra Tests
  --
  -------------------------------------------------------------------
  
  
  main :: IO ()
  main = defaultMain benchList
  
  benchList :: [Benchmark]
  benchList = fmap buildBench props
  
  props :: [(String,(Cl3 -> Bool))]
  props = [("Testing log.exp Identity:", prop_LogExp),
           ("Testing exp.log Identity:", prop_ExpLog),
           ("Testing abs*signum law:", prop_AbsSignum),
           ("Testing the definition of recip:", prop_RecipDef),
           ("Testing recip.recip Identity:", prop_RecipID),
           ("Testing sin.asin Identity:", prop_SinAsin),
           ("Testing asin.sin Identity:", prop_AsinSin),
           ("Testing cos.acos Identity:", prop_CosAcos),
           ("Testing acos.cos Identity:", prop_AcosCos),
           ("Testing sinh.asinh Identity:", prop_SinhAsinh),
           ("Testing asinh.sinh Identity:", prop_AsinhSinh),
           ("Testing cosh.acosh Identity:", prop_CoshAcosh),
           ("Testing acosh.cosh Identity:", prop_AcoshCosh),
           ("Testing acosh.cosh Identity2:", prop_AcoshCosh2),
           ("Testing Double Sin Identity:", prop_DubSin),
           ("Testing Double Cos Identity:", prop_DubCos),
           ("Testing Double Tan Identity:", prop_DubTan),
           ("Testing Double Sinh Identity:", prop_DubSinh),
           ("Testing Double Cosh Identity:", prop_DubCosh),
           ("Testing Double Tanh Identity:", prop_DubTanh),
           ("Testing Positive Sin Shift Identity:", prop_PosSinShift),
           ("Testing Negative Sin Shift Identity:", prop_NegSinShift),
           ("Testing sin^2+cos^2 Identity:", prop_SinSqCosSq),
           ("Testing cosh^2-sinh^2 Identity:", prop_CoshSqmSinhSq),
           ("Testing Symmetry of Cosh:", prop_SymCosh),
           ("Testing Symmetry of Sinh:", prop_SymSinh),
           ("Testing Double I Sin:", prop_DoubleISin)]
  
  buildBench :: (String,(Cl3 -> Bool)) -> Benchmark
  buildBench (name, prop) = runWithEnv $ \cliffs -> bench name (nfIO $ test cliffs)
    where
      test :: [Cl3] -> IO ()
      test ([]) = return ()
      test (cl:cls) =
        if prop cl
        then test cls
        else error $ "Failed on input: " ++ show cl
  
  runWithEnv :: ([Cl3] -> Benchmark) -> Benchmark
  runWithEnv = (env listRandCliffs)
  
  listRandCliffs :: IO [Cl3]
  listRandCliffs = do
    randCliff <-(replicateM 5000000).randomRIO $ (R 0, R 3)
    return (inputs ++ randCliff)
  
  -- Standard inputs and special cases of projectors and nilpotents
  inputs :: [Cl3]
  inputs = [R 0
           ,APS 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
           ,PV 0.5 0 0 0.5
           ,PV 0.5 0 0 (-0.5)
           ,BPV 0.5 0 0 0 (-0.5) 0
           ,BPV 0.5 0 0 0 0.5 0
           ,R 1
           ,R (-1)
           ,R pi
           ,R (pi/2)
           ,R (pi/4)
           ,V3 1 0 0
           ,APS 1 0.5 0 0 0 0.5 0 0
           ,APS 1 0.5 0 0 0 (-0.5) 0 0
           ,PV 1 1 0 0
           ,V3 1 0 0
           ,V3 (-1) 0 0
           ,V3 0 1 0
           ,V3 0 (-1) 0
           ,V3 0 0 1
           ,V3 0 0 (-1)
           ,V3 pi 0 0
           ,V3 (pi/2) 0 0
           ,V3 (pi/4) 0 0
           ,BV 1 0 0
           ,BV (-1) 0 0
           ,BV 0 1 0
           ,BV 0 (-1) 0
           ,BV 0 0 1
           ,BV 0 0 (-1)
           ,BV pi 0 0
           ,BV (pi/2) 0 0
           ,BV (pi/4) 0 0
           ,I 1
           ,I (-1)
           ,I pi
           ,I (pi/2)
           ,I (pi/4)
           ]
  
  
  -------------------------------------------------------
  -- | A set of properties to test
  -------------------------------------------------------
  
  prop_LogExp :: Cl3 -> Bool
  prop_LogExp (cliffor) = (abs cliffor > 10) || (
    let cliffor' = unWrapIPartEigs cliffor  -- imaginary part of log.exp repeats
  -- round off errors get large for exp larger than 5 use spectproj (log.exp) for accuracy
  -- note: +/- i*pi are not really poles but cause issues due to cancelation for (BV pi 0 0)
    in poles [I (-pi), I (pi)] cliffor' || (log (exp cliffor') ≈≈ cliffor'))
  
  -- log 0 is -Inf, Infinite vectors don't play nice
  -- spectproj (exp.log) doesn't have this issue
  prop_ExpLog :: Cl3 -> Bool
  prop_ExpLog (cliffor) = (lsv cliffor < tol) || (exp (log cliffor) ≈≈ cliffor)
  
  prop_AbsSignum :: Cl3 -> Bool
  prop_AbsSignum (cliffor) = abs cliffor * signum cliffor ≈≈ cliffor
  
  prop_RecipDef :: Cl3 -> Bool
  prop_RecipDef (cliffor) = (lsv cliffor < tol) || (recip cliffor * cliffor ≈≈ 1)
  
  -- singular inputs don't recip also suffers from roundoff errors at large values
  prop_RecipID :: Cl3 -> Bool
  prop_RecipID (cliffor) = (lsv cliffor < tol) || (recip (recip cliffor) ≈≈ cliffor)
  
  prop_SinAsin :: Cl3 -> Bool
  prop_SinAsin (cliffor) = if hasNilpotent cliffor
                           then poles [R 1, R (-1)] cliffor || (sin (asin cliffor) ≈≈ cliffor)
                           else sin (asin cliffor) ≈≈ cliffor
  
  prop_AsinSin :: Cl3 -> Bool
  prop_AsinSin (cliffor) = (abs cliffor > 10) || (asin (sin cliffor) ≈≈ (I (-1) * log (0.5 * (exp (I 1 * cliffor) - exp (mIx cliffor)) +
                                                                                       sqrt (1+0.25*(exp (mIx cliffor) - exp (I 1 * cliffor))^2))))
  
  prop_CosAcos :: Cl3 -> Bool
  prop_CosAcos (cliffor) = if hasNilpotent cliffor
                               then poles [R 1, R (-1)] cliffor || (cos (acos cliffor) ≈≈ cliffor)
                               else cos (acos cliffor) ≈≈ cliffor
  
  prop_AcosCos :: Cl3 -> Bool
  prop_AcosCos (cliffor) = (abs cliffor > 10) || (if hasNilpotent cliffor
                                                  then poles [R 0, pi, negate pi] cliffor || (acos (cos cliffor) ≈≈ 0.5 * (pi - 2 * asin(cos cliffor)))
                                                  else acos (cos cliffor) ≈≈ 0.5 * (pi - 2 * asin(cos cliffor)))
  
  prop_SinhAsinh :: Cl3 -> Bool
  prop_SinhAsinh (cliffor) = sinh (asinh cliffor) ≈≈ cliffor
  
  prop_AsinhSinh :: Cl3 -> Bool
  prop_AsinhSinh (cliffor) = (abs cliffor > 10) || (asinh (sinh cliffor) ≈≈ log (0.5*(exp cliffor - exp (negate cliffor)) +
                                                                                     sqrt (0.25 * (exp cliffor - exp (negate cliffor))^2 + 1)))
  
  prop_CoshAcosh :: Cl3 -> Bool
  prop_CoshAcosh (cliffor) = if hasNilpotent cliffor
                             then poles [R 1, R (-1)] cliffor || (cosh (acosh cliffor) ≈≈ cliffor)
                             else cosh (acosh cliffor) ≈≈ cliffor
  
  prop_AcoshCosh :: Cl3 -> Bool
  prop_AcoshCosh (cliffor) = acosh (cosh cliffor) ≈≈ log (0.5*(exp cliffor + exp (negate cliffor)) +
                                                          sqrt (0.5*(exp cliffor + exp (negate cliffor)) - 1) *
                                                          sqrt (0.5*(exp cliffor + exp (negate cliffor)) + 1))
  
  prop_AcoshCosh2 :: Cl3 -> Bool
  prop_AcoshCosh2 (cliffor) = acosh (cosh cliffor) ≈≈ log (cosh cliffor + sqrt (cosh cliffor - 1) * sqrt (cosh cliffor + 1))
  
  prop_DubSin :: Cl3 -> Bool
  prop_DubSin (cliffor) = sin (2 * cliffor) ≈≈ 2 * sin cliffor * cos cliffor
  
  prop_DubCos :: Cl3 -> Bool
  prop_DubCos (cliffor) = cos (2 * cliffor) ≈≈ cos cliffor ^ 2 - sin cliffor ^ 2
  
  prop_DubTan :: Cl3 -> Bool
  prop_DubTan (cliffor) = poles [R (-pi), R (-3*pi/4), R (-pi/2), R (-pi/4), R (pi/4), R (pi/2), R (3*pi/4), R (pi)] cliffor ||
                          (tan (2 * cliffor) ≈≈ (2 * tan cliffor) / (1 - tan cliffor ^ 2))
  
  prop_DubSinh :: Cl3 -> Bool
  prop_DubSinh (cliffor) = sinh (2 * cliffor) ≈≈ 2 * sinh cliffor * cosh cliffor
  
  prop_DubCosh :: Cl3 -> Bool
  prop_DubCosh (cliffor) = cosh (2 * cliffor) ≈≈ 2 * cosh cliffor ^ 2 - 1
  
  -- The test has poles at imaginary eigenvalues of n*pi/4 even is poles in the denominator and odd is poles in the numerator
  -- The poles are a source of a loss of precision.
  prop_DubTanh :: Cl3 -> Bool
  prop_DubTanh (cliffor) = poles [I (-pi), I (-3*pi/4), I (-pi/2), I (-pi/4), I (pi/4), I (pi/2), I (3*pi/4), I (pi)] cliffor ||
                           (tanh (2 * cliffor) ≈≈ (2 * tanh cliffor) / (1 + tanh cliffor ^ 2))
  
  prop_PosSinShift :: Cl3 -> Bool
  prop_PosSinShift (cliffor) = sin (pi/2 + cliffor) ≈≈ cos cliffor
  
  prop_NegSinShift :: Cl3 -> Bool
  prop_NegSinShift (cliffor) = sin (pi/2 - cliffor) ≈≈ cos cliffor
  
  prop_SinSqCosSq :: Cl3 -> Bool
  prop_SinSqCosSq (cliffor) = (abs cliffor > 10) || (sin cliffor ^ 2 + cos cliffor ^ 2 ≈≈ 1)
  
  prop_CoshSqmSinhSq :: Cl3 -> Bool
  prop_CoshSqmSinhSq (cliffor) = (abs cliffor > 10) || (cosh cliffor ^ 2 - sinh cliffor ^ 2 ≈≈ 1)
  
  prop_SymCosh :: Cl3 -> Bool
  prop_SymCosh (cliffor) = cosh (negate cliffor) ≈≈ cosh cliffor
  
  prop_SymSinh :: Cl3 -> Bool
  prop_SymSinh (cliffor) = sinh (negate cliffor) ≈≈ negate (sinh cliffor)
  
  prop_DoubleISin :: Cl3 -> Bool
  prop_DoubleISin (cliffor) = 2 * I 1 * sin cliffor ≈≈ exp(I 1 * cliffor) - exp (mIx cliffor)
  
  -- | Composition Sub-Algebras have a distributive norm over multiplication,
  -- like this:
  -- 
  -- > norm $ clif * clif' = norm clif * norm clif'
  --
  -- Strangly the constructor combinations with the "= True" don't play nice
  -- with 'abs' they are the constructors with non-zero zero-divisors.
  prop_CompAlg :: (Cl3, Cl3) -> Bool
  prop_CompAlg (PV{}, PV{}) = True
  prop_CompAlg (PV{}, BPV{}) = True
  prop_CompAlg (PV{}, TPV{}) = True
  prop_CompAlg (PV{}, APS{}) = True
  prop_CompAlg (BPV{}, PV{}) = True
  prop_CompAlg (TPV{}, PV{}) = True
  prop_CompAlg (APS{}, PV{}) = True
  prop_CompAlg (BPV{}, BPV{}) = True
  prop_CompAlg (BPV{}, TPV{}) = True
  prop_CompAlg (BPV{}, APS{}) = True
  prop_CompAlg (TPV{}, BPV{}) = True
  prop_CompAlg (APS{}, BPV{}) = True
  prop_CompAlg (TPV{}, TPV{}) = True
  prop_CompAlg (TPV{}, APS{}) = True
  prop_CompAlg (APS{}, TPV{}) = True
  prop_CompAlg (APS{}, APS{}) = True
  prop_CompAlg (cliffor, cliffor') = abs ( cliffor * cliffor') ≈≈ abs cliffor * abs cliffor'
  
  ----------------------------------------------------
  -- Helper functions for the properties
  ----------------------------------------------------
  
  -- | '≈≈' aproximately equal, using a mean squared error like calculation
  -- across the 8 dimensional vector space of APS.  The properties are 
  -- equivelent symbolicly but differ due to numerical errors.
  (≈≈) :: Cl3 -> Cl3 -> Bool
  (toAPS -> (APS a0 a1 a2 a3 a23 a31 a12 a123)) ≈≈ (toAPS -> (APS b0 b1 b2 b3 b23 b31 b12 b123)) =
    let m0 = (a0 - b0)^2
        m1 = (a1 - b1)^2
        m2 = (a2 - b2)^2
        m3 = (a3 - b3)^2
        m23 = (a23 - b23)^2
        m31 = (a31 - b31)^2
        m12 = (a12 - b12)^2
        m123 = (a123 - b123)^2
        sumsq = m0 + m1 + m2 + m3 + m23 + m31 + m12 + m123
        var = sumsq / 8
    in var <= 2e-13
  _ ≈≈ _ = error "Everything passed to (≈≈) should be caught by toAPS/APS pattern match"
  infix 4 ≈≈
  
  -- | 'poles' a function that tests if a cliffor is one of the defined poles
  poles :: [Cl3] -> Cl3 -> Bool
  poles [] _ = False
  poles [p] cliffor = eig1 `closeTo` p || eig2 `closeTo` p
    where (eig1,eig2) = eigvals cliffor
  poles (p:ps) cliffor = (eig1 `closeTo` p || eig2 `closeTo` p) || poles ps cliffor
    where (eig1,eig2) = eigvals cliffor
  
  -- | 'closeTo' used with poles to determine if an eigenvalue is close to a pole
  -- the current threshold is 1e-3
  closeTo :: Cl3 -> Cl3 -> Bool
  closeTo (toC -> (C a0 a123)) (toC -> (C b0 b123)) =
    let diffR = abs (a0 - b0)
        diffI = abs (a123 - b123)
        magDiff = sqrt (diffR^2 + diffI^2)
    in magDiff < 2e-3
  closeTo _ _ = error "Everything passed to 'closeTo' should be caught by toC/C pattern match"
  
  -- | 'unWrapIPartEigs' a function to reduce the magnitude of the imaginary
  -- portion of the Eigenvalues
  unWrapIPartEigs :: Cl3 -> Cl3
  unWrapIPartEigs cliffor = reduce $ spectraldcmp unWrapI id cliffor
    where unWrapI (R a0) = R a0
          unWrapI (I a123) | a123 > pi = unWrapI $ I (a123 - 2*pi)
                           | a123 < (-pi) = unWrapI $ I (a123 + 2*pi)
                           | otherwise = I a123
          unWrapI (C a0 a123) | a123 > pi = unWrapI $ C a0 (a123 - 2*pi)
                              | a123 < (-pi) = unWrapI $ C a0 (a123 + 2*pi)
                              | otherwise = C a0 a123
          unWrapI _ = error "unWrapI should only be unWrapping R I and C"
  
  -- End of File
  