{-# LANGUAGE UndecidableInstances #-}
--------------------------------------------------------------------------------
-- |
-- Module      :  HGeometry.Number.Real.Symbolic
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--
-- Defines a type 'Symbolic i r' that represents "indexed" numbers of
-- type r that have been "symbolically" perturbed by some arbitraryily
-- small \(\varepsilon\).
--
-- With indexed numbers we mean that every number has some index, and
-- the size of the pertubation \(\varepsilon\) is proportional to this index.
--
-- This is useful for implementing a more general form of "Simulation of Simplicity"
-- as described in:
--
--
--
--------------------------------------------------------------------------------
module HGeometry.Number.Real.Symbolic(
    EpsFold
  , eps, mkEpsFold
  , evalEps

  , hasNoPertubation
  , factors
  , suitableBase

  , Term(..), term, constantFactor

  , Symbolic
  , constant, symbolic, perturb

  , toTerms
  , signOf
  , roundToConstant

  , SoSRational, sosRational
  ) where

import           Control.Lens
import           Data.Foldable (toList)
import qualified Data.List as List
import qualified Data.Map.Merge.Strict as Map
import qualified Data.Map.Strict as Map
import           Data.Maybe (isNothing)
import           HGeometry.Number.Ratio.Generalized (GRatio, (%))
import           HGeometry.Sign (Sign(..))
-- 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 (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
$c<> :: forall i. Ord i => EpsFold i -> EpsFold i -> EpsFold i
<> :: EpsFold i -> EpsFold i -> EpsFold i
$csconcat :: forall i. Ord i => NonEmpty (EpsFold i) -> EpsFold i
sconcat :: NonEmpty (EpsFold i) -> EpsFold i
$cstimes :: forall i b. (Ord i, Integral b) => b -> EpsFold i -> EpsFold i
stimes :: forall b. Integral b => b -> 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
$cmempty :: forall i. Ord i => EpsFold i
mempty :: EpsFold i
$cmappend :: forall i. Ord i => EpsFold i -> EpsFold i -> EpsFold i
mappend :: EpsFold i -> EpsFold i -> EpsFold i
$cmconcat :: forall i. Ord i => [EpsFold i] -> EpsFold i
mconcat :: [EpsFold i] -> EpsFold i
Monoid)

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

-- | Creates the term \(\varepsilon(i)\)
eps :: i -> EpsFold i
eps :: forall i. 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

-- | Constructs an epsfold from a list of elements
mkEpsFold :: Ord i => [i] -> EpsFold i
mkEpsFold :: forall i. Ord i => [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 m a. Monoid m => (a -> m) -> [a] -> m
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 :: forall i. 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


-- | Evaluate an eps-fold using the choice of eps from the lemma above, i.e.
-- \(\varepsilon(i) = \delta^{d^i}\)
-- >>> evalEps 2 (1/2) $ mkEpsFold [1]
-- 0.25
-- >>> evalEps 2 (1/2) $ mkEpsFold [2]
-- 6.25e-2
-- >>> evalEps 2 (1/2) $ mkEpsFold [1,2]
-- 1.5625e-2
evalEps            :: (Fractional r, Integral i, Integral j) => j -> r -> EpsFold i -> r
evalEps :: forall r i j.
(Fractional r, Integral i, Integral j) =>
j -> r -> EpsFold i -> r
evalEps j
d r
delta EpsFold i
ef = r
delta r -> j -> r
forall a b. (Num a, Integral b) => a -> b -> a
^ [j] -> j
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ j
d j -> i -> j
forall a b. (Num a, Integral b) => a -> b -> a
^ i
i | i
i <- Bag i -> [i]
forall a. Bag a -> [a]
forall (t :: * -> *) a. Foldable t => t a -> [a]
toList (Bag i -> [i]) -> Bag i -> [i]
forall a b. (a -> b) -> a -> b
$ EpsFold i -> Bag i
forall i. EpsFold i -> Bag i
factors EpsFold i
ef]


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 a. Bag a -> [a]
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.


-- | Test if the epsfold has no pertubation at all (i.e. if it is \(\Pi_{\emptyset}\)
hasNoPertubation        :: EpsFold i -> Bool
hasNoPertubation :: forall i. EpsFold i -> Bool
hasNoPertubation (Pi Bag i
b) = Bag i -> Bool
forall a. Bag a -> 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
$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
/= :: Term i r -> Term i r -> Bool
Eq,(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
$cfmap :: forall i a b. (a -> b) -> Term i a -> Term i b
fmap :: forall a b. (a -> b) -> Term i a -> Term i b
$c<$ :: forall i a b. a -> Term i b -> Term i a
<$ :: forall a b. a -> Term i b -> Term i a
Functor)

-- | Lens to access the constant 'c' in the term.
constantFactor :: Lens' (Term i r) r
constantFactor :: forall i r (f :: * -> *).
Functor f =>
(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 :: forall r i. 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 Ordering -> Ordering -> Ordering
forall a. Semigroup a => a -> a -> a
<> r
c r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` r
d
                         -- note that the e2 and e1 are flipped;
                         -- i.e. if the numbers are negative this
                         -- reverses the ordering of the eps-values
                         -- involved.  (the c and d in the secondary
                         -- comparison are not flipped since they
                         -- should still incorporate the signs.)
                       (r
0,r
0)   -> Ordering
EQ -- both terms are zero.
                       (r
1,r
1)   -> EpsFold i
e1 EpsFold i -> EpsFold i -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` EpsFold i
e2 Ordering -> Ordering -> Ordering
forall a. Semigroup a => a -> a -> a
<> r
c r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` r
d

                       (-1,r
_)  -> Ordering
LT
                       (r
_,-1)  -> Ordering
GT

                       (r
0,r
1)   -> Ordering
LT
                       (r
1,r
_)   -> 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.


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

-- | Represents a Sum of terms, i.e. given a sequence of terms Term
-- \(c_i\) \(I_i\) s where \(c_i\) is a constant and \(I_i\) is a set
-- of indices, a Symbolic represents an expression of the form:
--
-- \[
--   \sum c_i \Pi_{j in I_i} \varepsilon(j)
-- \]
--
-- 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 ((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
$cfmap :: forall i a b. (a -> b) -> Symbolic i a -> Symbolic i b
fmap :: forall a b. (a -> b) -> Symbolic i a -> Symbolic i b
$c<$ :: forall i a b. a -> Symbolic i b -> Symbolic i a
<$ :: forall a b. a -> Symbolic i b -> Symbolic i a
Functor)


-- | Produces a list of terms, in decreasing order of significance
toTerms         :: Symbolic i r -> [Term i r]
toTerms :: forall i r. 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 :: forall r i. (Num r, Eq r) => 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

-- | Drops all the epsilon terms
roundToConstant :: Num r => Symbolic i r -> r
roundToConstant :: forall r i. Num r => Symbolic i r -> r
roundToConstant Symbolic i r
s = case Symbolic i r -> [Term i r]
forall i r. Symbolic i r -> [Term i r]
toTerms Symbolic i r
s of
                      (Term r
c EpsFold i
e:[Term i r]
_) | EpsFold i -> Bool
forall i. EpsFold i -> Bool
hasNoPertubation EpsFold i
e -> r
c
                      [Term i r]
_                                 -> r
0


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 a b. (a -> b) -> Symbolic i a -> Symbolic i b
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





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

-- | Creates a constant symbolic value
constant   :: Ord i => r -> Symbolic i r
constant :: forall i r. Ord i => 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 :: forall i r. Ord i => 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 :: forall r i. (Num r, Ord i) => 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 [ (i -> EpsFold i
forall i. i -> EpsFold i
eps i
i,r
1) , (EpsFold i
forall a. Monoid a => a
mempty,r
c) ]
  -- note that the empty epsfold is more significant than any other epsfold

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

-- | 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
$cshowsPrec :: forall a. Show a => Int -> Bag a -> ShowS
showsPrec :: Int -> Bag a -> ShowS
$cshow :: forall a. Show a => Bag a -> String
show :: Bag a -> String
$cshowList :: forall a. Show a => [Bag a] -> ShowS
showList :: [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
$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
/= :: 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
$ccompare :: forall a. Ord a => Bag a -> Bag a -> Ordering
compare :: Bag a -> Bag a -> Ordering
$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
>= :: Bag a -> Bag a -> Bool
$cmax :: forall a. Ord a => Bag a -> Bag a -> Bag a
max :: Bag a -> Bag a -> Bag a
$cmin :: forall a. Ord a => Bag a -> Bag a -> Bag a
min :: Bag a -> Bag a -> Bag a
Ord)

-- | Construct a singleton bag
singleton   :: k -> Bag k
singleton :: forall k. 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 :: forall m a. Monoid m => (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 m a. Monoid m => (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 :: forall a. 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 :: forall a. Ord a => 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


-- | Get the maximum in the bag
maximum'         :: Bag b -> Maybe b
maximum' :: forall b. Bag b -> Maybe b
maximum' (Bag Map b Int
m) = ((b, Int) -> b) -> Maybe (b, Int) -> Maybe b
forall a b. (a -> b) -> Maybe a -> 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 :: forall a. Bag a -> Int
maxMultiplicity (Bag Map a Int
m) = [Int] -> Int
forall a. Ord a => [a] -> a
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


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

-- | Number type representing
type SoSRational i r = GRatio (Symbolic i r)

-- | Helper to construct sosRationals
sosRational :: (Ord i, Eq r, Num r) => Symbolic i r -> Symbolic i r -> SoSRational i r
sosRational :: forall i r.
(Ord i, Eq r, Num r) =>
Symbolic i r -> Symbolic i r -> SoSRational i r
sosRational = Symbolic i r -> Symbolic i r -> GRatio (Symbolic i r)
forall a. (Eq a, Num a) => a -> a -> GRatio a
(%)

-- test :: GRatio Double
-- test = 15


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


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