--------------------------------------------------------------------------------
-- |
-- Module      :  Algorithms.Geometry.SoS.Symbolic
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--
--------------------------------------------------------------------------------
module Algorithms.Geometry.SoS.Symbolic(
    EpsFold
  , eps, mkEpsFold
  , hasNoPertubation
  , factors
  , suitableBase

  , Term(..), term, constantFactor

  , Symbolic
  , constant, symbolic, perturb

  , toTerms
  , signOf
  ) where

import           Algorithms.Geometry.SoS.Sign (Sign(..))
import           Control.Lens
import           Data.Foldable (toList)
import qualified Data.List as List
import qualified Data.Map as Map
import qualified Data.Map.Merge.Strict as Map
import           Data.Maybe (isNothing)
import           Test.QuickCheck (Arbitrary(..), listOf)
import           Test.QuickCheck.Instances ()

--------------------------------------------------------------------------------
-- * EpsFolds

{-
Let \(\mathcal{I}\) be a bag with indices, let \(c\) be an upper
bound on the number of times a single item may occur in
\(\mathcal{I}\), and let \(\varepsilon\) be a function mapping indices
to real numbers that satisfies:

1. \(0 < \varepsilon(j) < 1\), for all \(1 \leq j\),
2. \(\prod_{0 \leq i \leq j} \varepsilon(i)^c > \varepsilon(k)\), for all \(1 \leq j < k\)

Note that such a function exists:

\begin{lemma}
  \label{lem:condition_2}
  Let \(\delta \in (0,1)\) and \(d \geq c+1\). The function
  \(\varepsilon(i) = \delta^{d^i}\) satisfies condition 2.
\end{lemma}

\begin{proof}
  By transitivity it suffices to argue this for \(k=j+1\):

  \begin{align*}
           &\qquad \prod_{0 \leq i \leq j} \varepsilon(i)^c > \varepsilon(j+1) \\
    \equiv &\qquad \prod_{0 \leq i \leq j} (\delta^{d^i})^c > \delta^{d^{j+1}}\\
    \equiv &\qquad \prod_{0 \leq i \leq j} \delta^{cd^i}    > \delta^{d^{j+1}} \\
    \equiv &\qquad \delta^{\sum_{0 \leq i \leq j} cd^i} > \delta^{d^{j+1}} &
                                                                    \text{using
                                                                    }
                                                                    \delta \in (0,1)\\
    \equiv &\qquad \sum_{0 \leq i \leq j} cd^i < d^{j+1} \\
    \equiv &\qquad c\sum_{0 \leq i \leq j} d^i < d^{j+1} \\
  \end{align*}

  We prove this by induction.

  For the base case \(j=0\): we have \(0 < 1\), which is trivially true.

  For the step case we have the induction hypothesis
  \(c\sum_{0 \leq i \leq j} d^i < d^{j+1}\), and we have to prove that
  \(c\sum_{0 \leq i \leq j+1} d^i < d^{j+2}\):

  \begin{align*}
    c\sum_{0 \leq i \leq j+1} d^i
    &= cd^{j+1} + c\sum_{0 \leq i \leq j} d^i \\
    &< cd^{j+1} + d^{j+1}   & \text{using IH}  \\
    &= (c+1)d^{j+1}        & \text{using that } c+1 \leq d \\
    &\leq dd^{j+1}  \\
    &=d^{j+2}
  \end{align*}
  This completes the proof.
\end{proof}






An EpsFold now represents the term

\[ \prod_{i \in \mathcal{I}} \varepsilon(i) \]

for some bag \(\mathcal{I}\).


Let \(\mathcal{J}\) be some sub-bag of \(\mathcal{I}\). Note that
condition 2 implies that:

\(\prod_{i \in \mathcal{J}} \varepsilon(i) > \varepsilon(k)\), for all \(1 \leq j < k\)

This means that when comparing two EpsFolds, say \(e_1\) and \(e_2\),
representing bags \(\mathcal{I}_1\) and \(\mathcal{I}_2\),
respectively. It suffices to compare the largest index
\(j \in \mathcal{I}_1\setminus\mathcal{I}_2\) with the largest index
\(k \in \mathcal{I}_2\setminus\mathcal{I}_1\). We have that

\(e_1 > e_2\) if and only if \(j < k\).
-}
newtype EpsFold i = Pi (Bag i) deriving (b -> EpsFold i -> EpsFold i
NonEmpty (EpsFold i) -> EpsFold i
EpsFold i -> EpsFold i -> EpsFold i
(EpsFold i -> EpsFold i -> EpsFold i)
-> (NonEmpty (EpsFold i) -> EpsFold i)
-> (forall b. Integral b => b -> EpsFold i -> EpsFold i)
-> Semigroup (EpsFold i)
forall b. Integral b => b -> EpsFold i -> EpsFold i
forall i. Ord i => NonEmpty (EpsFold i) -> EpsFold i
forall i. Ord i => EpsFold i -> EpsFold i -> EpsFold i
forall i b. (Ord i, Integral b) => b -> EpsFold i -> EpsFold i
forall a.
(a -> a -> a)
-> (NonEmpty a -> a)
-> (forall b. Integral b => b -> a -> a)
-> Semigroup a
stimes :: b -> EpsFold i -> EpsFold i
$cstimes :: forall i b. (Ord i, Integral b) => b -> EpsFold i -> EpsFold i
sconcat :: NonEmpty (EpsFold i) -> EpsFold i
$csconcat :: forall i. Ord i => NonEmpty (EpsFold i) -> EpsFold i
<> :: EpsFold i -> EpsFold i -> EpsFold i
$c<> :: forall i. Ord i => EpsFold i -> EpsFold i -> EpsFold i
Semigroup,Semigroup (EpsFold i)
EpsFold i
Semigroup (EpsFold i)
-> EpsFold i
-> (EpsFold i -> EpsFold i -> EpsFold i)
-> ([EpsFold i] -> EpsFold i)
-> Monoid (EpsFold i)
[EpsFold i] -> EpsFold i
EpsFold i -> EpsFold i -> EpsFold i
forall i. Ord i => Semigroup (EpsFold i)
forall i. Ord i => EpsFold i
forall i. Ord i => [EpsFold i] -> EpsFold i
forall i. Ord i => EpsFold i -> EpsFold i -> EpsFold i
forall a.
Semigroup a -> a -> (a -> a -> a) -> ([a] -> a) -> Monoid a
mconcat :: [EpsFold i] -> EpsFold i
$cmconcat :: forall i. Ord i => [EpsFold i] -> EpsFold i
mappend :: EpsFold i -> EpsFold i -> EpsFold i
$cmappend :: forall i. Ord i => EpsFold i -> EpsFold i -> EpsFold i
mempty :: EpsFold i
$cmempty :: forall i. Ord i => EpsFold i
$cp1Monoid :: forall i. Ord i => Semigroup (EpsFold i)
Monoid)

-- | Gets the factors
factors         :: EpsFold i -> Bag i
factors :: EpsFold i -> Bag i
factors (Pi Bag i
is) = Bag i
is

-- | Creates the term \(\varepsilon(i)\)
eps :: i -> EpsFold i
eps :: i -> EpsFold i
eps = Bag i -> EpsFold i
forall i. Bag i -> EpsFold i
Pi (Bag i -> EpsFold i) -> (i -> Bag i) -> i -> EpsFold i
forall b c a. (b -> c) -> (a -> b) -> a -> c
. i -> Bag i
forall k. k -> Bag k
singleton

mkEpsFold :: Ord i => [i] -> EpsFold i
mkEpsFold :: [i] -> EpsFold i
mkEpsFold = Bag i -> EpsFold i
forall i. Bag i -> EpsFold i
Pi (Bag i -> EpsFold i) -> ([i] -> Bag i) -> [i] -> EpsFold i
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (i -> Bag i) -> [i] -> Bag i
forall (t :: * -> *) m a.
(Foldable t, Monoid m) =>
(a -> m) -> t a -> m
foldMap i -> Bag i
forall k. k -> Bag k
singleton



-- | computes a base 'd' that can be used as:
--
-- \( \varepsilon(i) = \varepsilon^{d^i} \)
suitableBase :: EpsFold i -> Int
suitableBase :: EpsFold i -> Int
suitableBase = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
2 (Int -> Int) -> (EpsFold i -> Int) -> EpsFold i -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
+) (Int -> Int) -> (EpsFold i -> Int) -> EpsFold i -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Bag i -> Int
forall a. Bag a -> Int
maxMultiplicity (Bag i -> Int) -> (EpsFold i -> Bag i) -> EpsFold i -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. EpsFold i -> Bag i
forall i. EpsFold i -> Bag i
factors

instance Show i => Show (EpsFold i) where
  showsPrec :: Int -> EpsFold i -> ShowS
showsPrec Int
d (Pi Bag i
b) = Bool -> ShowS -> ShowS
showParen (Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
app_prec) (ShowS -> ShowS) -> ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$
                         String -> ShowS
showString String
"Pi " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [i] -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
d (Bag i -> [i]
forall (t :: * -> *) a. Foldable t => t a -> [a]
toList Bag i
b)
    where
      app_prec :: Int
app_prec = Int
10


instance Ord i => Eq (EpsFold i) where
  EpsFold i
e1 == :: EpsFold i -> EpsFold i -> Bool
== EpsFold i
e2 = (EpsFold i
e1 EpsFold i -> EpsFold i -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` EpsFold i
e2) Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
EQ

instance Ord i => Ord (EpsFold i) where
  (Pi Bag i
e1) compare :: EpsFold i -> EpsFold i -> Ordering
`compare` (Pi Bag i
e2) = Maybe i
k Maybe i -> Maybe i -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` Maybe i
j -- note that k and j are flipped here
    where
      j :: Maybe i
j = Bag i -> Maybe i
forall b. Bag b -> Maybe b
maximum' (Bag i -> Maybe i) -> Bag i -> Maybe i
forall a b. (a -> b) -> a -> b
$ Bag i
e1 Bag i -> Bag i -> Bag i
forall a. Ord a => Bag a -> Bag a -> Bag a
`difference` Bag i
e2
      k :: Maybe i
k = Bag i -> Maybe i
forall b. Bag b -> Maybe b
maximum' (Bag i -> Maybe i) -> Bag i -> Maybe i
forall a b. (a -> b) -> a -> b
$ Bag i
e2 Bag i -> Bag i -> Bag i
forall a. Ord a => Bag a -> Bag a -> Bag a
`difference` Bag i
e1
    -- note: If the terms are all the same, the difference of the bags is empty
    -- and thus both e1e2 and e2e1 are Nothing and thus equal.

    -- otherwise, let j be the largest term that is in e1 but not in e2.
    -- If e2 does not have any terms at all (Nothing) it will be bigger than e1
    --
    -- if e2 does have a term, let k be the largest one, then the
    -- biggest of those terms is the pair whose indices comes first.

instance (Arbitrary i, Ord i) => Arbitrary (EpsFold i) where
  arbitrary :: Gen (EpsFold i)
arbitrary = [i] -> EpsFold i
forall i. Ord i => [i] -> EpsFold i
mkEpsFold ([i] -> EpsFold i) -> ([i] -> [i]) -> [i] -> EpsFold i
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [i] -> [i]
forall a. Int -> [a] -> [a]
take Int
4 ([i] -> EpsFold i) -> Gen [i] -> Gen (EpsFold i)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Gen i -> Gen [i]
forall a. Gen a -> Gen [a]
listOf Gen i
forall a. Arbitrary a => Gen a
arbitrary


-- | Test if the epsfold has no pertubation at all (i.e. if it is \(\Pi_{\emptyset}\)
hasNoPertubation        :: EpsFold i -> Bool
hasNoPertubation :: EpsFold i -> Bool
hasNoPertubation (Pi Bag i
b) = Bag i -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null Bag i
b


--------------------------------------------------------------------------------
-- * Terms

-- | A term 'Term c es' represents a term:
--
-- \[ c \Pi_{i \in es} \varepsilon(i)
-- \]
--
-- for a constant c and an arbitrarily small value \(\varepsilon\),
-- parameterized by i.
data Term i r = Term r (EpsFold i) deriving (Term i r -> Term i r -> Bool
(Term i r -> Term i r -> Bool)
-> (Term i r -> Term i r -> Bool) -> Eq (Term i r)
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
forall i r. (Ord i, Eq r) => Term i r -> Term i r -> Bool
/= :: Term i r -> Term i r -> Bool
$c/= :: forall i r. (Ord i, Eq r) => Term i r -> Term i r -> Bool
== :: Term i r -> Term i r -> Bool
$c== :: forall i r. (Ord i, Eq r) => Term i r -> Term i r -> Bool
Eq,a -> Term i b -> Term i a
(a -> b) -> Term i a -> Term i b
(forall a b. (a -> b) -> Term i a -> Term i b)
-> (forall a b. a -> Term i b -> Term i a) -> Functor (Term i)
forall a b. a -> Term i b -> Term i a
forall a b. (a -> b) -> Term i a -> Term i b
forall i a b. a -> Term i b -> Term i a
forall i a b. (a -> b) -> Term i a -> Term i b
forall (f :: * -> *).
(forall a b. (a -> b) -> f a -> f b)
-> (forall a b. a -> f b -> f a) -> Functor f
<$ :: a -> Term i b -> Term i a
$c<$ :: forall i a b. a -> Term i b -> Term i a
fmap :: (a -> b) -> Term i a -> Term i b
$cfmap :: forall i a b. (a -> b) -> Term i a -> Term i b
Functor)

-- | Lens to access the constant 'c' in the term.
constantFactor :: Lens' (Term i r) r
constantFactor :: (r -> f r) -> Term i r -> f (Term i r)
constantFactor = (Term i r -> r)
-> (Term i r -> r -> Term i r) -> Lens (Term i r) (Term i r) r r
forall s a b t. (s -> a) -> (s -> b -> t) -> Lens s t a b
lens (\(Term r
c EpsFold i
_) -> r
c) (\(Term r
_ EpsFold i
es) r
c -> r -> EpsFold i -> Term i r
forall i r. r -> EpsFold i -> Term i r
Term r
c EpsFold i
es)


instance (Show i, Show r) => Show (Term i r) where
  showsPrec :: Int -> Term i r -> ShowS
showsPrec Int
d (Term r
c EpsFold i
es) = Bool -> ShowS -> ShowS
showParen (Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
up_prec) (ShowS -> ShowS) -> ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$
                               Int -> r -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec (Int
up_prec Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) r
c
                             ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> ShowS
showString String
" * "
                             ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> EpsFold i -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec (Int
up_prec Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) EpsFold i
es
    where
      up_prec :: Int
up_prec = Int
5


-- | Creates a singleton term
term     :: r -> i -> Term i r
term :: r -> i -> Term i r
term r
r i
i = r -> EpsFold i -> Term i r
forall i r. r -> EpsFold i -> Term i r
Term r
r (EpsFold i -> Term i r) -> EpsFold i -> Term i r
forall a b. (a -> b) -> a -> b
$ i -> EpsFold i
forall i. i -> EpsFold i
eps i
i

instance (Ord i, Ord r, Num r) => Ord (Term i r) where
  (Term r
c EpsFold i
e1) compare :: Term i r -> Term i r -> Ordering
`compare` (Term r
d EpsFold i
e2) = case (EpsFold i -> Bool
forall i. EpsFold i -> Bool
hasNoPertubation EpsFold i
e1, EpsFold i -> Bool
forall i. EpsFold i -> Bool
hasNoPertubation EpsFold i
e2) of
                                        (Bool
True,Bool
True) -> r
c    r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` r
d
                                        (Bool, Bool)
_           -> case (r -> r
forall a. Num a => a -> a
signum r
c, r -> r
forall a. Num a => a -> a
signum r
d) of
                                                         (-1,-1) -> EpsFold i
e2 EpsFold i -> EpsFold i -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` EpsFold i
e1
                                                         (r
0,r
0)   -> EpsFold i
e1 EpsFold i -> EpsFold i -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` EpsFold i
e2
                                                         (r
1,r
1)   -> EpsFold i
e1 EpsFold i -> EpsFold i -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` EpsFold i
e2
                                                         (-1,r
_)  -> Ordering
LT
                                                         (r
_,-1)  -> Ordering
GT
                                                         (r, r)
_       -> String -> Ordering
forall a. HasCallStack => String -> a
error String
"SoS: Term.ord absurd"
  -- If both the eps folds are zero, and thus we just have constants
  -- then we should compare the individual terms.

  -- if *one* of the two has an eps term, then we can choose eps to be
  -- arbitrarily small, i.e. small enough so that that terms is
  -- actually smaller than the other term.  this is reflected since
  -- findMax will then return a Noting, which is smaller than anything
  -- else

  -- if both terms have epsilon terms, we first look at the sign. If
  -- they have non-negative signs we compare the eps-folds as in the
  -- paper. (Lemma 3.3). If both are negative, that reverses the
  -- ordering. If the signs are different then we can base the
  -- ordering on that.

instance (Arbitrary r, Arbitrary (EpsFold i), Ord i) => Arbitrary (Term i r) where
  arbitrary :: Gen (Term i r)
arbitrary = r -> EpsFold i -> Term i r
forall i r. r -> EpsFold i -> Term i r
Term (r -> EpsFold i -> Term i r)
-> Gen r -> Gen (EpsFold i -> Term i r)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Gen r
forall a. Arbitrary a => Gen a
arbitrary Gen (EpsFold i -> Term i r) -> Gen (EpsFold i) -> Gen (Term i r)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Gen (EpsFold i)
forall a. Arbitrary a => Gen a
arbitrary

--------------------------------------------------------------------------------
-- * Symbolic

-- | Represents a Sum of terms, i.e. a value that has the form:
--
-- \[
--   \sum c \Pi_i \varepsilon(i)
-- \]
--
-- The terms are represented in order of decreasing significance.
--
-- The main idea in this type is that, if symbolic values contains
-- \(\varepsilon(i)\) terms we can always order them. That is, two
-- Symbolic terms will be equal only if:
--
-- - they contain *only* a constant term (that is equal)
-- - they contain the exact same \(\varepsilon\)-fold.
--
newtype Symbolic i r = Sum (Map.Map (EpsFold i) r) deriving (a -> Symbolic i b -> Symbolic i a
(a -> b) -> Symbolic i a -> Symbolic i b
(forall a b. (a -> b) -> Symbolic i a -> Symbolic i b)
-> (forall a b. a -> Symbolic i b -> Symbolic i a)
-> Functor (Symbolic i)
forall a b. a -> Symbolic i b -> Symbolic i a
forall a b. (a -> b) -> Symbolic i a -> Symbolic i b
forall i a b. a -> Symbolic i b -> Symbolic i a
forall i a b. (a -> b) -> Symbolic i a -> Symbolic i b
forall (f :: * -> *).
(forall a b. (a -> b) -> f a -> f b)
-> (forall a b. a -> f b -> f a) -> Functor f
<$ :: a -> Symbolic i b -> Symbolic i a
$c<$ :: forall i a b. a -> Symbolic i b -> Symbolic i a
fmap :: (a -> b) -> Symbolic i a -> Symbolic i b
$cfmap :: forall i a b. (a -> b) -> Symbolic i a -> Symbolic i b
Functor)

-- | Produces a list of terms, in decreasing order of significance
toTerms         :: Symbolic i r -> [Term i r]
toTerms :: Symbolic i r -> [Term i r]
toTerms (Sum Map (EpsFold i) r
m) = ((EpsFold i, r) -> Term i r) -> [(EpsFold i, r)] -> [Term i r]
forall a b. (a -> b) -> [a] -> [b]
map (\(EpsFold i
i,r
c) -> r -> EpsFold i -> Term i r
forall i r. r -> EpsFold i -> Term i r
Term r
c EpsFold i
i) ([(EpsFold i, r)] -> [Term i r])
-> (Map (EpsFold i) r -> [(EpsFold i, r)])
-> Map (EpsFold i) r
-> [Term i r]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Map (EpsFold i) r -> [(EpsFold i, r)]
forall k a. Map k a -> [(k, a)]
Map.toDescList (Map (EpsFold i) r -> [Term i r])
-> Map (EpsFold i) r -> [Term i r]
forall a b. (a -> b) -> a -> b
$ Map (EpsFold i) r
m

-- | Computing the Sign of an expression. (Nothing represents zero)
signOf   :: (Num r, Eq r) => Symbolic i r -> Maybe Sign
signOf :: Symbolic i r -> Maybe Sign
signOf Symbolic i r
e = case (r -> Bool) -> [r] -> [r]
forall a. (a -> Bool) -> [a] -> [a]
List.dropWhile (r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
0) ([r] -> [r]) -> ([Term i r] -> [r]) -> [Term i r] -> [r]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Term i r -> r) -> [Term i r] -> [r]
forall a b. (a -> b) -> [a] -> [b]
map (\(Term r
c EpsFold i
_) -> r -> r
forall a. Num a => a -> a
signum r
c) ([Term i r] -> [r]) -> [Term i r] -> [r]
forall a b. (a -> b) -> a -> b
$ Symbolic i r -> [Term i r]
forall i r. Symbolic i r -> [Term i r]
toTerms Symbolic i r
e of
             []     -> Maybe Sign
forall a. Maybe a
Nothing
             (-1:[r]
_) -> Sign -> Maybe Sign
forall a. a -> Maybe a
Just Sign
Negative
             [r]
_      -> Sign -> Maybe Sign
forall a. a -> Maybe a
Just Sign
Positive

instance (Ord i, Eq r, Num r) => Eq (Symbolic i r) where
  Symbolic i r
e1 == :: Symbolic i r -> Symbolic i r -> Bool
== Symbolic i r
e2 = Maybe Sign -> Bool
forall a. Maybe a -> Bool
isNothing (Maybe Sign -> Bool) -> Maybe Sign -> Bool
forall a b. (a -> b) -> a -> b
$ Symbolic i r -> Maybe Sign
forall r i. (Num r, Eq r) => Symbolic i r -> Maybe Sign
signOf (Symbolic i r
e1 Symbolic i r -> Symbolic i r -> Symbolic i r
forall a. Num a => a -> a -> a
- Symbolic i r
e2)

instance (Ord i, Ord r, Num r) => Ord (Symbolic i r) where
  Symbolic i r
e1 compare :: Symbolic i r -> Symbolic i r -> Ordering
`compare` Symbolic i r
e2 = case Symbolic i r -> Maybe Sign
forall r i. (Num r, Eq r) => Symbolic i r -> Maybe Sign
signOf (Symbolic i r
e1 Symbolic i r -> Symbolic i r -> Symbolic i r
forall a. Num a => a -> a -> a
- Symbolic i r
e2) of
                      Maybe Sign
Nothing       -> Ordering
EQ
                      Just Sign
Negative -> Ordering
LT
                      Just Sign
Positive -> Ordering
GT

instance (Ord i, Num r, Eq r) => Num (Symbolic i r) where
  (Sum Map (EpsFold i) r
e1) + :: Symbolic i r -> Symbolic i r -> Symbolic i r
+ (Sum Map (EpsFold i) r
e2) = Map (EpsFold i) r -> Symbolic i r
forall i r. Map (EpsFold i) r -> Symbolic i r
Sum (Map (EpsFold i) r -> Symbolic i r)
-> Map (EpsFold i) r -> Symbolic i r
forall a b. (a -> b) -> a -> b
$ SimpleWhenMissing (EpsFold i) r r
-> SimpleWhenMissing (EpsFold i) r r
-> SimpleWhenMatched (EpsFold i) r r r
-> Map (EpsFold i) r
-> Map (EpsFold i) r
-> Map (EpsFold i) r
forall k a c b.
Ord k =>
SimpleWhenMissing k a c
-> SimpleWhenMissing k b c
-> SimpleWhenMatched k a b c
-> Map k a
-> Map k b
-> Map k c
Map.merge SimpleWhenMissing (EpsFold i) r r
forall (f :: * -> *) k x. Applicative f => WhenMissing f k x x
Map.preserveMissing -- insert things only in e1
                                        SimpleWhenMissing (EpsFold i) r r
forall (f :: * -> *) k x. Applicative f => WhenMissing f k x x
Map.preserveMissing -- insert things only in e2
                                        SimpleWhenMatched (EpsFold i) r r r
forall k. WhenMatched Identity k r r r
combine
                                        Map (EpsFold i) r
e1 Map (EpsFold i) r
e2
    where
      -- if things are in both e1 and e2, we add the constant terms. If they are non-zero
      -- we use this value in the map. Otherwise we drop it.
      combine :: WhenMatched Identity k r r r
combine = (k -> r -> r -> Maybe r) -> WhenMatched Identity k r r r
forall (f :: * -> *) k x y z.
Applicative f =>
(k -> x -> y -> Maybe z) -> WhenMatched f k x y z
Map.zipWithMaybeMatched
                (\k
_ r
c r
d -> let x :: r
x = r
c r -> r -> r
forall a. Num a => a -> a -> a
+ r
d in if r
x r -> r -> Bool
forall a. Eq a => a -> a -> Bool
/= r
0 then r -> Maybe r
forall a. a -> Maybe a
Just r
x else Maybe r
forall a. Maybe a
Nothing)
    -- Symbolic $ Map.unionWith (+) ts ts'

  negate :: Symbolic i r -> Symbolic i r
negate = (r -> r) -> Symbolic i r -> Symbolic i r
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap r -> r
forall a. Num a => a -> a
negate

  (Sum Map (EpsFold i) r
ts) * :: Symbolic i r -> Symbolic i r -> Symbolic i r
* (Sum Map (EpsFold i) r
ts') = Map (EpsFold i) r -> Symbolic i r
forall i r. Map (EpsFold i) r -> Symbolic i r
Sum (Map (EpsFold i) r -> Symbolic i r)
-> Map (EpsFold i) r -> Symbolic i r
forall a b. (a -> b) -> a -> b
$ (r -> r -> r) -> [(EpsFold i, r)] -> Map (EpsFold i) r
forall k a. Ord k => (a -> a -> a) -> [(k, a)] -> Map k a
Map.fromListWith r -> r -> r
forall a. Num a => a -> a -> a
(+) [ (EpsFold i
es EpsFold i -> EpsFold i -> EpsFold i
forall a. Semigroup a => a -> a -> a
<> EpsFold i
es',r
cr -> r -> r
forall a. Num a => a -> a -> a
*r
d)
                                                    | (EpsFold i
es, r
c) <- Map (EpsFold i) r -> [(EpsFold i, r)]
forall k a. Map k a -> [(k, a)]
Map.toList Map (EpsFold i) r
ts
                                                    , (EpsFold i
es',r
d) <- Map (EpsFold i) r -> [(EpsFold i, r)]
forall k a. Map k a -> [(k, a)]
Map.toList Map (EpsFold i) r
ts'
                                                    , r
cr -> r -> r
forall a. Num a => a -> a -> a
*r
d r -> r -> Bool
forall a. Eq a => a -> a -> Bool
/= r
0
                                                    ]

  fromInteger :: Integer -> Symbolic i r
fromInteger Integer
x = r -> Symbolic i r
forall i r. Ord i => r -> Symbolic i r
constant (Integer -> r
forall a. Num a => Integer -> a
fromInteger Integer
x)

  signum :: Symbolic i r -> Symbolic i r
signum Symbolic i r
s = case Symbolic i r -> Maybe Sign
forall r i. (Num r, Eq r) => Symbolic i r -> Maybe Sign
signOf Symbolic i r
s of
               Maybe Sign
Nothing       -> Symbolic i r
0
               Just Sign
Negative -> (-Symbolic i r
1)
               Just Sign
Positive -> Symbolic i r
1

  abs :: Symbolic i r -> Symbolic i r
abs Symbolic i r
x | Symbolic i r -> Symbolic i r
forall a. Num a => a -> a
signum Symbolic i r
x Symbolic i r -> Symbolic i r -> Bool
forall a. Eq a => a -> a -> Bool
== -Symbolic i r
1 = (-Symbolic i r
1)Symbolic i r -> Symbolic i r -> Symbolic i r
forall a. Num a => a -> a -> a
*Symbolic i r
x
        | Bool
otherwise      = Symbolic i r
x


instance (Show i, Show r) => Show (Symbolic i r) where
  showsPrec :: Int -> Symbolic i r -> ShowS
showsPrec Int
d Symbolic i r
s = Bool -> ShowS -> ShowS
showParen (Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
app_prec) (ShowS -> ShowS) -> ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$
                    String -> ShowS
showString String
"Sum " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [Term i r] -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
d (Symbolic i r -> [Term i r]
forall i r. Symbolic i r -> [Term i r]
toTerms Symbolic i r
s)
    where
      app_prec :: Int
app_prec = Int
10

instance (Arbitrary r, Ord i, Arbitrary (EpsFold i)) => Arbitrary (Symbolic i r) where
  arbitrary :: Gen (Symbolic i r)
arbitrary = Map (EpsFold i) r -> Symbolic i r
forall i r. Map (EpsFold i) r -> Symbolic i r
Sum (Map (EpsFold i) r -> Symbolic i r)
-> Gen (Map (EpsFold i) r) -> Gen (Symbolic i r)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Gen (Map (EpsFold i) r)
forall a. Arbitrary a => Gen a
arbitrary

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

-- | Creates a constant symbolic value
constant   :: Ord i => r -> Symbolic i r
constant :: r -> Symbolic i r
constant r
c = Map (EpsFold i) r -> Symbolic i r
forall i r. Map (EpsFold i) r -> Symbolic i r
Sum (Map (EpsFold i) r -> Symbolic i r)
-> Map (EpsFold i) r -> Symbolic i r
forall a b. (a -> b) -> a -> b
$ EpsFold i -> r -> Map (EpsFold i) r
forall k a. k -> a -> Map k a
Map.singleton EpsFold i
forall a. Monoid a => a
mempty r
c

-- | Creates a symbolic vlaue with a single indexed term. If you just need a constant (i.e. non-indexed), use 'constant'
symbolic     :: Ord i => r -> i -> Symbolic i r
symbolic :: r -> i -> Symbolic i r
symbolic r
r i
i = Map (EpsFold i) r -> Symbolic i r
forall i r. Map (EpsFold i) r -> Symbolic i r
Sum (Map (EpsFold i) r -> Symbolic i r)
-> Map (EpsFold i) r -> Symbolic i r
forall a b. (a -> b) -> a -> b
$ EpsFold i -> r -> Map (EpsFold i) r
forall k a. k -> a -> Map k a
Map.singleton (i -> EpsFold i
forall i. i -> EpsFold i
eps i
i) r
r

-- | given the value c and the index i, creates the perturbed value
-- \(c + \varepsilon(i)\)
perturb      :: (Num r, Ord i) => r -> i -> Symbolic i r
perturb :: r -> i -> Symbolic i r
perturb r
c i
i = Map (EpsFold i) r -> Symbolic i r
forall i r. Map (EpsFold i) r -> Symbolic i r
Sum (Map (EpsFold i) r -> Symbolic i r)
-> Map (EpsFold i) r -> Symbolic i r
forall a b. (a -> b) -> a -> b
$ [(EpsFold i, r)] -> Map (EpsFold i) r
forall k a. Eq k => [(k, a)] -> Map k a
Map.fromAscList [ (EpsFold i
forall a. Monoid a => a
mempty,r
c) , (i -> EpsFold i
forall i. i -> EpsFold i
eps i
i,r
1) ]


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

-- | The word specifiies how many *duplicates* there are. I.e. If the
-- Bag maps k to i, then k has multiplicity i+1.
newtype Bag a = Bag (Map.Map a Int) deriving (Int -> Bag a -> ShowS
[Bag a] -> ShowS
Bag a -> String
(Int -> Bag a -> ShowS)
-> (Bag a -> String) -> ([Bag a] -> ShowS) -> Show (Bag a)
forall a. Show a => Int -> Bag a -> ShowS
forall a. Show a => [Bag a] -> ShowS
forall a. Show a => Bag a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Bag a] -> ShowS
$cshowList :: forall a. Show a => [Bag a] -> ShowS
show :: Bag a -> String
$cshow :: forall a. Show a => Bag a -> String
showsPrec :: Int -> Bag a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> Bag a -> ShowS
Show,Bag a -> Bag a -> Bool
(Bag a -> Bag a -> Bool) -> (Bag a -> Bag a -> Bool) -> Eq (Bag a)
forall a. Eq a => Bag a -> Bag a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Bag a -> Bag a -> Bool
$c/= :: forall a. Eq a => Bag a -> Bag a -> Bool
== :: Bag a -> Bag a -> Bool
$c== :: forall a. Eq a => Bag a -> Bag a -> Bool
Eq,Eq (Bag a)
Eq (Bag a)
-> (Bag a -> Bag a -> Ordering)
-> (Bag a -> Bag a -> Bool)
-> (Bag a -> Bag a -> Bool)
-> (Bag a -> Bag a -> Bool)
-> (Bag a -> Bag a -> Bool)
-> (Bag a -> Bag a -> Bag a)
-> (Bag a -> Bag a -> Bag a)
-> Ord (Bag a)
Bag a -> Bag a -> Bool
Bag a -> Bag a -> Ordering
Bag a -> Bag a -> Bag a
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
forall a. Ord a => Eq (Bag a)
forall a. Ord a => Bag a -> Bag a -> Bool
forall a. Ord a => Bag a -> Bag a -> Ordering
forall a. Ord a => Bag a -> Bag a -> Bag a
min :: Bag a -> Bag a -> Bag a
$cmin :: forall a. Ord a => Bag a -> Bag a -> Bag a
max :: Bag a -> Bag a -> Bag a
$cmax :: forall a. Ord a => Bag a -> Bag a -> Bag a
>= :: Bag a -> Bag a -> Bool
$c>= :: forall a. Ord a => Bag a -> Bag a -> Bool
> :: Bag a -> Bag a -> Bool
$c> :: forall a. Ord a => Bag a -> Bag a -> Bool
<= :: Bag a -> Bag a -> Bool
$c<= :: forall a. Ord a => Bag a -> Bag a -> Bool
< :: Bag a -> Bag a -> Bool
$c< :: forall a. Ord a => Bag a -> Bag a -> Bool
compare :: Bag a -> Bag a -> Ordering
$ccompare :: forall a. Ord a => Bag a -> Bag a -> Ordering
$cp1Ord :: forall a. Ord a => Eq (Bag a)
Ord,Gen (Bag a)
Gen (Bag a) -> (Bag a -> [Bag a]) -> Arbitrary (Bag a)
Bag a -> [Bag a]
forall a. (Ord a, Arbitrary a) => Gen (Bag a)
forall a. (Ord a, Arbitrary a) => Bag a -> [Bag a]
forall a. Gen a -> (a -> [a]) -> Arbitrary a
shrink :: Bag a -> [Bag a]
$cshrink :: forall a. (Ord a, Arbitrary a) => Bag a -> [Bag a]
arbitrary :: Gen (Bag a)
$carbitrary :: forall a. (Ord a, Arbitrary a) => Gen (Bag a)
Arbitrary)

singleton   :: k -> Bag k
singleton :: k -> Bag k
singleton k
x = Map k Int -> Bag k
forall a. Map a Int -> Bag a
Bag (Map k Int -> Bag k) -> Map k Int -> Bag k
forall a b. (a -> b) -> a -> b
$ k -> Int -> Map k Int
forall k a. k -> a -> Map k a
Map.singleton k
x Int
0


instance Foldable Bag where
  -- ^ Takes multiplicity into account.
  foldMap :: (a -> m) -> Bag a -> m
foldMap a -> m
f (Bag Map a Int
m) =
    (a -> Int -> m) -> Map a Int -> m
forall m k a. Monoid m => (k -> a -> m) -> Map k a -> m
Map.foldMapWithKey (\a
k Int
d -> (a -> m) -> [a] -> m
forall (t :: * -> *) m a.
(Foldable t, Monoid m) =>
(a -> m) -> t a -> m
foldMap a -> m
f (Int -> a -> [a]
forall a. Int -> a -> [a]
List.replicate (Int -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) a
k)) Map a Int
m
  null :: Bag a -> Bool
null (Bag Map a Int
m) = Map a Int -> Bool
forall k a. Map k a -> Bool
Map.null Map a Int
m

instance Ord k => Semigroup (Bag k) where
  (Bag Map k Int
m) <> :: Bag k -> Bag k -> Bag k
<> (Bag Map k Int
m') = Map k Int -> Bag k
forall a. Map a Int -> Bag a
Bag (Map k Int -> Bag k) -> Map k Int -> Bag k
forall a b. (a -> b) -> a -> b
$ (Int -> Int -> Int) -> Map k Int -> Map k Int -> Map k Int
forall k a. Ord k => (a -> a -> a) -> Map k a -> Map k a -> Map k a
Map.unionWith (\Int
d Int
d' -> Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d' Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Map k Int
m Map k Int
m'

instance Ord k => Monoid (Bag k) where
  mempty :: Bag k
mempty = Map k Int -> Bag k
forall a. Map a Int -> Bag a
Bag Map k Int
forall k a. Map k a
Map.empty

-- | Computes the difference of the two maps
difference                   :: Ord a => Bag a -> Bag a -> Bag a
difference :: Bag a -> Bag a -> Bag a
difference (Bag Map a Int
m1) (Bag Map a Int
m2) = Map a Int -> Bag a
forall a. Map a Int -> Bag a
Bag (Map a Int -> Bag a) -> Map a Int -> Bag a
forall a b. (a -> b) -> a -> b
$ (Int -> Int -> Maybe Int) -> Map a Int -> Map a Int -> Map a Int
forall k a b.
Ord k =>
(a -> b -> Maybe a) -> Map k a -> Map k b -> Map k a
Map.differenceWith Int -> Int -> Maybe Int
forall a. (Num a, Ord a) => a -> a -> Maybe a
updateCount Map a Int
m1 Map a Int
m2
  where
    updateCount :: a -> a -> Maybe a
updateCount a
i a
j = let d :: a
d = a
i a -> a -> a
forall a. Num a => a -> a -> a
- a
j -- note that we should actually compare (i+1) and (j+1)
                      in if a
d a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 then Maybe a
forall a. Maybe a
Nothing -- we have no copies left
                                   else a -> Maybe a
forall a. a -> Maybe a
Just (a -> Maybe a) -> a -> Maybe a
forall a b. (a -> b) -> a -> b
$ a
d a -> a -> a
forall a. Num a => a -> a -> a
- a
1


maximum'         :: Bag b -> Maybe b
maximum' :: Bag b -> Maybe b
maximum' (Bag Map b Int
m) = ((b, Int) -> b) -> Maybe (b, Int) -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (b, Int) -> b
forall a b. (a, b) -> a
fst (Maybe (b, Int) -> Maybe b)
-> (Map b Int -> Maybe (b, Int)) -> Map b Int -> Maybe b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Map b Int -> Maybe (b, Int)
forall k a. Map k a -> Maybe (k, a)
Map.lookupMax (Map b Int -> Maybe b) -> Map b Int -> Maybe b
forall a b. (a -> b) -> a -> b
$ Map b Int
m


-- | maximum multiplicity of an element in the bag
maxMultiplicity         :: Bag a -> Int
maxMultiplicity :: Bag a -> Int
maxMultiplicity (Bag Map a Int
m) = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ([Int] -> Int) -> (Map a Int -> [Int]) -> Map a Int -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int
0Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:) ([Int] -> [Int]) -> (Map a Int -> [Int]) -> Map a Int -> [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
+) ([Int] -> [Int]) -> (Map a Int -> [Int]) -> Map a Int -> [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Map a Int -> [Int]
forall k a. Map k a -> [a]
Map.elems (Map a Int -> Int) -> Map a Int -> Int
forall a b. (a -> b) -> a -> b
$ Map a Int
m