Approx.hs

{-|
    Module      :  Data.Number.ER.Real.Approx
    Description :  classes abstracting exact reals
    Copyright   :  (c) Michal Konecny
    License     :  BSD3

    Maintainer  :  mik@konecny.aow.cz
    Stability   :  experimental
    Portability :  portable

    Definitions of classes that describe what is
    required from arbitrary precision approximations
    of exact real numbers.
    
    We introduce two levels of abstraction for these
    approximations:
    
        * 'ERApprox' = 
            a *set* of approximated numbers whose size is
            measured using some fixed measure
        
        * 'ERIntApprox' = 
            an *interval* of real numbers with finitely
            representable endpoints 
    
    To be imported qualified, usually with the synonym RA.
-}
module Data.Number.ER.Real.Approx
(
    ERApprox(..),
    ERIntApprox(..),
    bounds2ira,
    effIx2ra,
    splitIRA,
--    checkShrinking,
    eqSingletons,
    leqSingletons,
    ltSingletons,
    equalIntervals,
    exactMiddle,
    maxExtensionR2R
)
where

import Data.Number.ER.BasicTypes
import qualified Data.Number.ER.ExtendedInteger as EI

import Data.Typeable

{-|
   A type whose elements represent sets that can be used
   to approximate a single extended real number with arbitrary precision.
-}
class (Fractional ra) => ERApprox ra where
    getPrecision :: ra -> Precision 
    {-^ 
            Precision is a measure of the set size.
            
            The default interpretation:
            
            * If the diameter of the set is d, then the precision
            should be near floor(- log_2 d).
    -}
    getGranularity :: ra -> Granularity
    -- ^ the lower the granularity the bigger the rounding errors
    setGranularity :: Granularity -> ra -> ra
    -- ^ increase or safely decrease granularity
    setMinGranularity :: Granularity -> ra -> ra
    -- ^ ensure granularity is not below the first arg
    isEmpty :: ra -> Bool 
    -- ^ true if this represents a computational error
    isBottom :: ra -> Bool 
    -- ^ true if this holds no information
    isExact :: ra -> Bool 
    -- ^ true if this is a singleton
    isDisjoint :: ra -> ra -> Bool
    isDisjoint a b = isEmpty $ a /\ b
    isInteriorDisjoint :: ra -> ra -> Bool
    isInteriorDisjoint a b = 
        isEmpty isect || isExact isect 
        where
        isect = a /\ b
    isBounded :: ra -> Bool 
    -- ^ true if the approximation excludes infinity
    bottomApprox :: ra 
    -- ^ the bottom element - any number
    emptyApprox :: ra 
    -- ^ the top element - error
    refines :: ra -> ra -> Bool 
    -- ^ first arg is a subset of the second arg
    (/\) :: ra -> ra -> ra 
    -- ^ join; combining two approximations of the same number
    intersectMeasureImprovement ::
        EffortIndex -> ra -> ra -> (ra, ra)
    {-^ 
            Like intersection but the second component:
            
             * measures improvement of the intersection relative to the first of the two approximations
             
             * is a positive number: 1 means no improvement, 2 means doubled precision, etc. 
    -}
    equalReals :: ra -> ra -> Maybe Bool 
    -- ^ nothing if overlapping and not singletons
    compareReals :: ra -> ra -> Maybe Ordering
    -- ^ nothing if overlapping and not singletons
    leqReals :: ra -> ra -> Maybe Bool
    -- ^ nothing if overlapping on interior or by a wrong endpoint
    equalApprox :: ra -> ra -> Bool
    -- ^ syntactic comparison
    compareApprox :: ra -> ra -> Ordering
    -- ^ syntactic linear ordering
    double2ra :: Double -> ra
    showApprox :: 
        Int {-^ number of relevant decimals to show -} ->
        Bool {-^ should show granularity -} ->
        Bool {-^ should show internal representation details -} ->
        ra {-^ the approximation to show -} ->
        String
    
{-|
    For a finite sequence of real approximations, determine
    whether it is a shrinking sequence.
-}    
checkShrinking ::
    (ERApprox ra) =>
    [ra] -> 
    Maybe (ra, ra)
checkShrinking [] = Nothing
checkShrinking [_] = Nothing
checkShrinking (a : b : rest) 
    | b `refines` a = checkShrinking (b : rest)
    | otherwise = Just (a,b)

        
{-|
   A type whose elements represent sets that can be used
   to approximate a recursive set of closed extended real number intervals 
   with arbitrary precision.
-}
--class (ERApprox sra) => SetOfRealsApprox sra where
--    (\/) :: sra -> sra -> sra -- ^ union; either approximation could be correct

{-|
   A type whose elements represent real *intervals* that can be used
   to approximate a single extended real number with arbitrary precision.

   Sometimes, these types can be used to approximate 
   a closed extended real number interval with arbitrary precision.
   Nevetheless, this is not guaranteed.
-}
class (ERApprox ira) => ERIntApprox ira 
    where
    doubleBounds :: ira -> (Double, Double) 
    floatBounds :: ira -> (Float, Float)
    integerBounds :: ira -> (EI.ExtendedInteger, EI.ExtendedInteger)
    bisectDomain :: 
        Maybe ira {-^ point to split at -} -> 
        ira {-^ interval to split -} -> 
        (ira, ira) -- ^ left and right, overlapping on a singleton
    defaultBisectPt :: ira -> ira
    -- | returns thin approximations of endpoints, in natural order 
    bounds :: ira -> (ira, ira)
    {-|
         meet, usually constructing interval from approximations of its endpoints
         
         This does not need to be the meet of the real intervals 
         but it has to be a maximal element in the set of all
         ira elements that are below the two parameters.
    -}
    (\/) :: ira -> ira -> ira
    
{-| 
    Inverse of 'bounds'.
-}
bounds2ira ::
    (ERIntApprox ira) =>
    (ira, ira) -> ira
bounds2ira (a,b) = a \/ b
    
{-|
    Assuming the arguments are singletons, equality is decidable.
-}
eqSingletons :: (ERApprox ra) => ra -> ra -> Bool
eqSingletons s1 s2 =  
    case equalReals s1 s2 of 
        Just b -> b
        _ -> False 

{-|
    Assuming the arguments are singletons, @<=@ is decidable.
-}
leqSingletons :: (ERApprox ra) => ra -> ra -> Bool
leqSingletons s1 s2 =  
    case compareReals s1 s2 of 
        Just EQ -> True
        Just LT -> True
        _ -> False 
        
{-|
    Assuming the arguments are singletons, @<@ is decidable.
-}
ltSingletons :: (ERApprox ra) => ra -> ra -> Bool
ltSingletons s1 s2 =  
    case compareReals s1 s2 of 
        Just LT -> True
        _ -> False 
        
{-|
    Return true if and only if the two intervals have equal endpoints.
-}
equalIntervals ::
    (ERIntApprox ira) => ira -> ira -> Bool
equalIntervals d1 d2 =
    d1L == d2L && d1U == d2U
    where
    (==) = eqSingletons
    (d1L, d1U) = bounds d1
    (d2L, d2U) = bounds d2


{-|    
    This function converts
    an effort index to a real number approximation.
    
    Useful when an effort index is used in a formula
    mixed with real approximations.  
-}
effIx2ra :: 
    (ERApprox ra) =>
    EffortIndex -> ra
effIx2ra = fromInteger . toInteger

{-|
    Split an interval to a sequence of intervals whose union is the
    original interval using a given sequence of cut points.
    The cut points are expected to be in increasing order and contained
    in the given interval.  Violations of this rule are tolerated.
-}
splitIRA ::
    (ERIntApprox ira) =>
    ira {-^ an interval to be split -} -> 
    [ira] {-^ approximations of the cut points in increasing order -} -> 
    [ira]
splitIRA interval splitPoints =
    doSplit [] end pointsRev
    where
    (start, end) = bounds interval
    pointsRev = reverse $ start : splitPoints
    doSplit previousSegments nextRight [] = previousSegments
    doSplit previousSegments nextRight (nextLeft : otherPoints) =
        doSplit (nextLeft \/ nextRight : previousSegments) nextLeft otherPoints

{-|
    * Return the endpoints of the interval as well as the exact midpoint.
    
    * To be able to do this, there may be a need to increase granularity.
    
    * All three singleton intervals are set to the same new granularity.
-}        
exactMiddle ::
    (ERIntApprox ira) =>
    ira ->
    (ira,ira,ira,Granularity)
exactMiddle dom =
    case isExact domM of
        True ->
            (domL, domM, domR, gran)
        False ->
            (domLhg, domMhg, domRhg, higherGran)
    where
    (domL, domR) = bounds dom
    gran = max (getGranularity domL) (getGranularity domR)
    domM = (domL + domR) / 2
    higherGran = gran + 1
    domLhg = setMinGranularity higherGran domL
    domRhg = setMinGranularity higherGran domR
    domMhg = (domLhg + domRhg) / 2
     
        
{-| 
    This produces a function that computes the maximal extension of the
    given function.  A maximal extension function has the property:
    f(I) = { f(x) | x in I }.  Here we get this property only for the
    limit function for its 'EffortIndex' tending to infinity.
-}
maxExtensionR2R ::
    (ERIntApprox ira) =>
    (EffortIndex -> ira -> [ira]) 
        {-^ returns a safe approximation of all extrema within the interval -} ->
    (EffortIndex -> ira -> ira) 
        {-^ a function behaving well on sequences that intersect to a point -} ->
    (EffortIndex -> ira -> ira) 
        {- ^ a function behaving well on sequences that intersect to a non-empty interval -}
maxExtensionR2R getExtremes f ix x
    | getPrecision x < effIx2prec ix =
        (f ix xL) \/ (f ix xR) \/ 
        (foldl (\/) emptyApprox $ getExtremes ix x)
    -- x is thin enough (?), don't bother evaluating by endpoints and extrema:
    | otherwise =
        f ix x
    where
    (xL, xR) = bounds x