Interval.hs

{-# OPTIONS_GHC -fno-warn-missing-methods #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE UndecidableInstances   #-}
{-# LANGUAGE FlexibleInstances   #-}
{-# LANGUAGE DeriveDataTypeable   #-}

{-|
    Module      :  Data.Number.ER.RnToRm.UnitDom.Approx.Interval
    Description :  arbitrary precision function enclosures on @[-1,1]^n@
    Copyright   :  (c) Michal Konecny
    License     :  BSD3

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

    A construction of an enclosure of a real function on
    the domain [-1,1]^n for some n using elements of some
    base (eg rational functions or polynomials).
-}
module Data.Number.ER.RnToRm.UnitDom.Approx.Interval 
(
    ERFnInterval(..),
    ERFnContext(..)
)
where

import qualified Data.Number.ER.Real.Base as B
import Data.Number.ER.Real.Approx.Interval
import Data.Number.ER.Real.Arithmetic.Elementary

import qualified Data.Number.ER.RnToRm.Approx as FA
import qualified Data.Number.ER.RnToRm.UnitDom.Approx as UFA
import qualified Data.Number.ER.RnToRm.UnitDom.Base as UFB
import qualified Data.Number.ER.Real.Approx as RA
import qualified Data.Number.ER.Real.Approx.Elementary as RAEL

import qualified Data.Number.ER.Real.DomainBox as DBox
import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)
import Data.Number.ER.BasicTypes

import Data.Number.ER.Misc

import qualified Data.Map as Map

import Data.Typeable
import Data.Generics.Basics
import Data.Binary

{- only for testing in ghci, to be removed: -}
--import Data.Number.ER.Real.DefaultRepr
--import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom
--import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.PolynomBase
--type FAPU = ERFnInterval (ERChebPoly B) IRA
--fapuConst1 = (UFA.const 0 [1]) :: FAPU
--fapuConst2 = (UFA.const 0 [2]) :: FAPU
{- end of testing specific code -}

data ERFnInterval fb ra =
    ERFnIntervalAny 
    {
        erfnContext :: ERFnContext
    }
    |
    ERFnInterval 
    {
        erfnUpper :: fb,
        erfnLowerNeg :: fb,
        erfnContext :: ERFnContext,
        erfnGlobal :: ra
    }
    deriving (Typeable, Data)

instance (Binary a, Binary b) => Binary (ERFnInterval a b) where
  put (ERFnIntervalAny a) = putWord8 0 >> put a
  put (ERFnInterval a b c d) = putWord8 1 >> put a >> put b >> put c >> put d
  get = do
    tag_ <- getWord8
    case tag_ of
      0 -> get >>= \a -> return (ERFnIntervalAny a)
      1 -> get >>= \a -> get >>= \b -> get >>= \c -> get >>= \d -> return (ERFnInterval a b c d)
      _ -> fail "no parse"
    

data ERFnContext =
    ERFnContext
    {
        erfnMaxDegree :: Int,
        erfnCoeffGranularity :: Granularity
    }
    deriving (Show, Typeable, Data)
    
instance Binary ERFnContext where
  put (ERFnContext a b) = put a >> put b
  get = get >>= \a -> get >>= \b -> return (ERFnContext a b)
    
    
erfnContextDefault =
    ERFnContext
    {
        erfnMaxDegree = 2,
        erfnCoeffGranularity = 10
    }
    
erfnContextUnify (ERFnContext dg1 gr1) (ERFnContext dg2 gr2) =
    ERFnContext (max dg1 dg2) (max gr1 gr2)

    
instance 
    (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
    Show (ERFnInterval fb ra)
    where
    show (ERFnIntervalAny _) = "ERFnIntervalAny"
    show (ERFnInterval h ln ctxt gl) =
        "\nERFnInterval"
        ++ "\n  upper = " ++ show h
        ++ "\n  lower = " ++ show (-ln)
--        ++ "  global = " ++ show gl ++ "\n"
--        ++ "  context = " ++ show ctxt ++ "\n"

instance
    (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
    Eq (ERFnInterval fb ra)
    where
    (ERFnInterval h1 ln1 ctxt1 gl1) 
            == (ERFnInterval h2 ln2 ctxt2 gl2) =
        error "ERFnInterval: equality not implemented"
    _ == _ =
        error "ERFnInterval: equality not implemented"

instance 
    (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
    Ord (ERFnInterval fb ra) 
    where
    compare 
            (ERFnInterval h1 ln1 ctxt1 gl1) 
            (ERFnInterval h2 ln2 ctxt2 gl2) =
        error "ERFnInterval: comparison not implemented; consider leqReals from class ERApprox instead"
    compare _ _ =
        error "ERFnInterval: comparison not implemented; consider leqReals from class ERApprox instead"
    
    
instance 
    (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
    Num (ERFnInterval fb ra)
    where
    fromInteger n = UFA.const [fromInteger n]
    negate f@(ERFnIntervalAny _) = f
    negate (ERFnInterval h ln ctxt gl) =
        (ERFnInterval ln h ctxt (negate gl))
    (ERFnInterval h1 ln1 ctxt1 gl1) + (ERFnInterval h2 ln2 ctxt2 gl2) =
        ERFnInterval (h1 + h2) (ln1 + ln2) ctxt (gl1 + gl2)
        where
        ctxt = erfnContextUnify ctxt1 ctxt2
    f1 + f2 = ERFnIntervalAny ctxt
        where
        ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)
    (ERFnInterval h1 ln1 ctxt1 gl1) * (ERFnInterval h2 ln2 ctxt2 gl2) =
        ERFnInterval h ln ctxt (gl1 * gl2)
        where
        (h, ln) =
            case (RA.leqReals 0 gl1, RA.leqReals gl1 0, RA.leqReals 0 gl2, RA.leqReals gl2 0) of
                (Just True, _, Just True, _) -> -- both non-negative
                    (h1h2, l1l2Neg)
                (_, Just True, _, Just True) -> -- both non-positive
                    (l1l2, h1h2Neg)
                (Just True, _, _, Just True) -> -- first non-negative, second non-positive
                    (l1h2, h1l2Neg)
                (_, Just True, Just True, _) -> -- first non-positive, second non-negative
                    (h1l2, l1h2Neg)
                _ -> -- one of both may be crossing zero
                    ((h1h2 `maxP` l1l2) `maxP` (h1l2 `maxP` l1h2),
                     (h1h2Neg `maxP` l1l2Neg) `maxP` (h1l2Neg `maxP` l1h2Neg))
                where
                h1h2 = UFB.reduceDegreeUp maxDegr $ h1 * h2
                h1h2Neg = UFB.reduceDegreeUp maxDegr $ (negate h1) * h2
                l1l2 = UFB.reduceDegreeUp maxDegr $ ln1 * ln2
                l1l2Neg = UFB.reduceDegreeUp maxDegr $ (negate ln1) * ln2
                h1l2 = UFB.reduceDegreeUp maxDegr $ h1 * (negate ln2)
                h1l2Neg = UFB.reduceDegreeUp maxDegr $ h1 * ln2
                l1h2 = UFB.reduceDegreeUp maxDegr $ (negate ln1) * h2
                l1h2Neg = UFB.reduceDegreeUp maxDegr $ ln1 * h2
                maxP p1 p2 = fst $ UFB.max maxDegr p1 p2
                     
        ctxt = erfnContextUnify ctxt1 ctxt2
        maxDegr = erfnMaxDegree ctxt
    f1 * f2 = ERFnIntervalAny ctxt
        where
        ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)
        
instance 
    (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
    Fractional (ERFnInterval fb ra) 
    where
    fromRational r = UFA.const [fromRational r]
    recip f@(ERFnIntervalAny _) = f
    recip (ERFnInterval h ln ctxt gl) 
        | certainNoZero =
            ERFnInterval lRecipUp hnRecipUp ctxt (recip gl)
        | otherwise = ERFnIntervalAny ctxt
        where
        certainNoZero =
            certainAboveZero || certainBelowZero
        certainAboveZero =
             UFB.upperBound ix ln < 0
        certainBelowZero =         
             UFB.upperBound ix h < 0 
        hnRecipUp =
            UFB.recipUp maxDegr ix (negate h) 
        lRecipUp =
            UFB.recipUp maxDegr ix (negate ln)
        maxDegr = erfnMaxDegree ctxt
        ix = int2effIx $ 3 * maxDegr         

instance
    (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
    RA.ERApprox (ERFnInterval fb ra) 
    where
    getGranularity (ERFnIntervalAny ctxt) = erfnCoeffGranularity ctxt
    getGranularity (ERFnInterval h ln ctxt gl) =
        max (erfnCoeffGranularity ctxt) $ 
            max (UFB.getGranularity h) (UFB.getGranularity ln)
    setGranularity gran (ERFnIntervalAny ctxt) = 
        ERFnIntervalAny $ ctxt { erfnCoeffGranularity = gran }
    setGranularity gran (ERFnInterval h ln ctxt gl) =
        ERFnInterval 
            (UFB.setGranularity gran h) (UFB.setGranularity gran ln) 
            (ctxt { erfnCoeffGranularity = gran }) gl
    setMinGranularity gran (ERFnIntervalAny ctxt) = 
        ERFnIntervalAny
            (ctxt { erfnCoeffGranularity = max gran (erfnCoeffGranularity ctxt) })
    setMinGranularity gran (ERFnInterval h ln ctxt gl) =
        ERFnInterval 
            (UFB.setMinGranularity gran h) (UFB.setMinGranularity gran ln) 
            (ctxt { erfnCoeffGranularity = max gran (erfnCoeffGranularity ctxt) }) gl
--    getPrecision (ERFnIntervalAny _) = 0
--    getPrecision f = intLog 2 (1 + (fst $ RA.integerBounds (FA.volume f))) -- wrong! 
    (ERFnInterval h1 ln1 ctxt1 gl1) /\ (ERFnInterval h2 ln2 ctxt2 gl2) =
        ERFnInterval (snd $ UFB.min maxDegr h1 h2) (snd $ UFB.min maxDegr ln1 ln2) ctxt (gl1 RA./\ gl2)
        where
        ctxt = erfnContextUnify ctxt1 ctxt2
        maxDegr = erfnMaxDegree ctxt
    (ERFnIntervalAny ctxt1) /\ (ERFnInterval h2 ln2 ctxt2 gl2) =
        ERFnInterval h2 ln2 ctxt gl2
        where
        ctxt = erfnContextUnify ctxt1 ctxt2
    (ERFnInterval h1 ln1 ctxt1 gl1) /\ (ERFnIntervalAny ctxt2) =
        ERFnInterval h1 ln1 ctxt gl1
        where
        ctxt = erfnContextUnify ctxt1 ctxt2
    f1 /\ f2 = ERFnIntervalAny ctxt
        where
        ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)
    leqReals = erfnintLeq

erfnintLeq left right
    | left `isClearlyBelow` right = Just True
    | right `isClearlyStrictlyBelow` left = Just False
    | otherwise = Nothing
    where
    isClearlyBelow (ERFnIntervalAny _) _ = False
    isClearlyBelow _ (ERFnIntervalAny _) = False
    isClearlyBelow f g
        | UFB.upperBound 10 (erfnUpper f + erfnLowerNeg g) <= 0 = True
        | otherwise = False
    isClearlyStrictlyBelow (ERFnIntervalAny _) _ = False
    isClearlyStrictlyBelow _ (ERFnIntervalAny _) = False
    isClearlyStrictlyBelow f g
        | UFB.upperBound 10 (erfnUpper f + erfnLowerNeg g) < 0 = True
        | otherwise = False

instance
    (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
    RA.ERIntApprox (ERFnInterval fb ra) 
    where
--    doubleBounds = :: ira -> (Double, Double) 
--    floatBounds :: ira -> (Float, Float)
--    integerBounds :: ira -> (ExtendedInteger, ExtendedInteger)
    bisectDomain maybePt (ERFnIntervalAny c) =
        error "ERFnInterval: RA.bisectDomain: cannot bisect ERFnIntervalAny"
    bisectDomain maybePt (ERFnInterval u ln c g) =
        (ERFnInterval midUp ln c g,
         ERFnInterval u (negate midDown) c g)
         where
         (midDown, midUp) =
            case maybePt of
                Nothing ->
                    (negate $ (ln - u) / 2, (u - ln) / 2)
                Just (ERFnInterval uPt lnPt _ _) ->
                    (negate lnPt, uPt)
    bounds (ERFnIntervalAny c) =
        error "ERFnInterval: RA.bounds: cannot get bounds for ERFnIntervalAny"
    bounds (ERFnInterval u ln c g) =
        (ERFnInterval (negate ln) ln c g,
         ERFnInterval u (negate u) c g) 
    (ERFnInterval u1 ln1 c1 g1) \/ (ERFnInterval u2 ln2 c2 g2) =
        ERFnInterval u ln c (g1 RA.\/ g2)
        where
        u = UFB.maxUp maxDegree u1 u2
        ln = UFB.maxUp maxDegree ln1 ln2
        c = erfnContextUnify c1 c2
        maxDegree = erfnMaxDegree c
    (ERFnIntervalAny ctxt1) \/ (ERFnInterval h2 ln2 ctxt2 gl2) =
        ERFnIntervalAny ctxt
        where
        ctxt = erfnContextUnify ctxt1 ctxt2
    (ERFnInterval h1 ln1 ctxt1 gl1) \/ (ERFnIntervalAny ctxt2) =
        ERFnIntervalAny ctxt
        where
        ctxt = erfnContextUnify ctxt1 ctxt2
    f1 \/ f2 = ERFnIntervalAny ctxt
        where
        ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)

instance
    (UFB.ERUnitFnBase boxb boxra varid b ra fb, RAEL.ERApproxElementary ra, RealFrac b) =>
    RAEL.ERApproxElementary (ERFnInterval fb ra) 
    where
    -- default abs does not work because we do not have Prelude.abs
    abs _ f@(ERFnIntervalAny _) = f
    abs _ (ERFnInterval u ln c g) =
        ERFnInterval maxulnUp maxunl0Dn c (abs g)
        where
        maxDegree = erfnMaxDegree c
        maxulnUp = snd $ UFB.max maxDegree u ln 
        maxunl0Dn =
            fst $ UFB.max maxDegree 0 $
                fst $ UFB.max maxDegree (- u) (- ln)
    exp ix f@(ERFnIntervalAny _) = f
    exp ix (ERFnInterval u ln c g) =
        ERFnInterval uExp lExpNeg c (RAEL.exp ix g)
        where
        maxDegree = erfnMaxDegree c
--        ix = int2effIx maxDegree
        uExp = snd $ UFB.exp maxDegree ix u
        lExpNeg = 
            negate $ fst $ UFB.exp maxDegree ix (negate ln) 
    sin ix f@(ERFnIntervalAny c) = 
        ERFnInterval 1 1 c ((-1) RA.\/ 1)
    sin ix (ERFnInterval u ln c g) =
--        unsafePrint
--        (
--            "ERFnInterval: RAEL.sin: "
--            ++ "\n u = " ++ show u
--            ++ "\n ln = " ++ show ln
--            ++ "\n uSin = " ++ show uSin
--            ++ "\n lSinNeg = " ++ show lSinNeg
--        ) $
        ERFnInterval uSin (- lSin) c (RAEL.sin ix g)
        where
        (lSin, uSin) = sincos True maxDegree ix u (-ln)  
        maxDegree = erfnMaxDegree c
    cos ix f@(ERFnIntervalAny c) =
        ERFnInterval 1 1 c ((-1) RA.\/ 1)
    cos ix (ERFnInterval u ln c g) =
--        unsafePrint
--        (
--            "ERFnInterval: RAEL.cos: "
--            ++ "\n u = " ++ show u
--            ++ "\n ln = " ++ show ln
--            ++ "\n uCos = " ++ show uCos
--            ++ "\n lCosNeg = " ++ show lCosNeg
--        ) $
        ERFnInterval uCos (- lCos) c (RAEL.cos ix g)
        where
        (lCos, uCos) = sincos False maxDegree ix u (-ln) 
        maxDegree = erfnMaxDegree c
    atan ix f@(ERFnIntervalAny c) =
        ERFnInterval 1 1 c ((-1) RA.\/ 1)
    atan ix (ERFnInterval u ln c g) =
--        unsafePrint
--        (
--            "ERFnInterval: RAEL.atan: "
--            ++ "\n u = " ++ show u
--            ++ "\n ln = " ++ show ln
--            ++ "\n uAtan = " ++ show uAtan
--            ++ "\n lAtanNeg = " ++ show lAtanNeg
--        ) $
        ERFnInterval uAtan lAtanNeg c (RAEL.atan ix g)
        where
        maxDegree = erfnMaxDegree c
--        ix = int2effIx maxDegree
        uAtan = snd $ UFB.atan maxDegree ix u
        lAtanNeg = 
            negate $ fst $ UFB.atan maxDegree ix (negate ln) 

sincos ::
    (UFB.ERUnitFnBase boxb boxra varid b ra fb, RAEL.ERApproxElementary ra, RealFrac b) =>
    Bool {-^ True iff sine, False iff cosine -} -> 
    Int {-^ maximum representation degree -} -> 
    EffortIndex {-^ how hard to try to eliminate truncation errors -} -> 
    fb ->
    fb ->
    (fb, fb)
sincos isSine maxDegree ix u l
    -- p - 2k*pi range within [-pi/2, pi/2]: 
    | ranfNear0 `RA.refines` plusMinusPiHalf =
--        unsafePrint
--        (
--            "ERFnInterval: sincos: [-pi/2, pi/2]: "
--            ++ "\n u = " ++ show u
--            ++ "\n l = " ++ show l
--            ++ "\n ranf = " ++ show ranf
--            ++ "\n k = " ++ show k
--            ++ "\n ranfNear0 = " ++ show ranfNear0
--        ) $
        case isSine of
            True -> sineShifted (- k2pi)
            False -> cosineShifted (- k2pi)
    -- p - 2k*pi range within [0, pi]: 
    | (ranfNear0 - piHalf) `RA.refines` plusMinusPiHalf =
--        unsafePrint
--        (
--            "ERFnInterval: sincos: [0, pi]: "
--            ++ "\n u = " ++ show u
--            ++ "\n l = " ++ show l
--            ++ "\n ranf = " ++ show ranf
--            ++ "\n k = " ++ show k
--            ++ "\n ranfNear0 = " ++ show ranfNear0
--        ) $
        case isSine of
            -- use sin(x) = cos(x - pi/2) and cos(x) = - sin(x - pi/2):
            True -> cosineShifted (- k2pi - piHalf)
            False -> sineShiftedNegated (- k2pi - piHalf)
    -- p - 2k*pi range within [-pi, 0]: 
    | (ranfNear0 + piHalf) `RA.refines` plusMinusPiHalf =
--        unsafePrint
--        (
--            "ERFnInterval: sincos: [-pi, 0]: "
--            ++ "\n u = " ++ show u
--            ++ "\n l = " ++ show l
--            ++ "\n ranf = " ++ show ranf
--            ++ "\n k = " ++ show k
--            ++ "\n ranfNear0 = " ++ show ranfNear0
--        ) $
        case isSine of
            -- use sin(x) = - cos(x + pi/2) and cos(x) = sin(x + pi/2):
            True -> cosineShiftedNegated (-k2pi + piHalf)
            False -> sineShifted (-k2pi + piHalf)
    -- p - 2k*pi range within [pi/2, 3pi/2]: 
    | (ranfNear0 - pi) `RA.refines` plusMinusPiHalf =
--        unsafePrint
--        (
--            "ERFnInterval: sincos: [pi/2, 3pi/2]: "
--            ++ "\n u = " ++ show u
--            ++ "\n l = " ++ show l
--            ++ "\n ranf = " ++ show ranf
--            ++ "\n k = " ++ show k
--            ++ "\n ranfNear0 = " ++ show ranfNear0
--        ) $
        -- use sin(x) = - sin(x - pi) and cos(x) = - cos(x - pi)
        case isSine of
            True -> sineShiftedNegated (- k2pi - pi)
            False -> cosineShiftedNegated (- k2pi - pi)
    | otherwise = 
--        unsafePrint
--        (
--            "ERFnInterval: sincos: big range: "
--            ++ "\n u = " ++ show u
--            ++ "\n l = " ++ show l
--            ++ "\n ranf = " ++ show ranf
--            ++ "\n k = " ++ show k
--            ++ "\n ranfNear0 = " ++ show ranfNear0
--        ) $
        (UFB.const (-1), UFB.const 1)
--    (expDownwards, expUpwards + valueAtRDnNeg + (UFB.const expRUp))
    where
    ranfNear0 = ranf - k2pi
    k2pi = k * 2 * pi
    plusMinusPiHalf = (-piHalfLO) RA.\/ piHalfLO
    pi = RAEL.pi ix  
    piHalf = pi / 2
    (piHalfLO, piHalfHI) = RA.bounds piHalf
    ranf = 
        ERInterval 
            (UFB.lowerBound 10 l) 
            (UFB.upperBound 10 u)
    k = 
        fromInteger $ floor $ 
            case (pi + ranf) / (2 * pi) of ERInterval lo hi -> lo

    sineShiftedNegated shift =
        boundsNegate $ sineShifted shift
        
    cosineShiftedNegated shift =
        boundsNegate $ cosineShifted shift

    boundsNegate (pLO, pHI) = (- pHI, - pLO)
        
    sineShifted shift =
        boundsAddErr shiftWidthB (lSinDown, uSinUp)
        where
        lSinDown = fst $ UFB.sin maxDegree ix (l `plusUp` shiftPoly)
        uSinUp = snd $ UFB.sin maxDegree ix (u `plusDown` shiftPoly)  
        shiftPoly = UFB.const shiftLOB
        ERInterval shiftLOB shiftHIB = shift
        shiftWidthB = shiftHIB - shiftLOB
    
    cosineShifted shift =
        boundsAddErr shiftWidthB $ 
            (UFB.minDown maxDegree lCosDown uCosDown,
             UFB.maxUp maxDegree lCosUp uCosUp 
                + (snd $ UFB.scale 0.5 (u-l))) -- important near 0
        where
        (lCosDown, lCosUp) = UFB.cos maxDegree ix (l `plusUp` shiftPoly)
        (uCosDown, uCosUp) = UFB.cos maxDegree ix (u `plusDown` shiftPoly)  
        shiftPoly = UFB.const shiftLOB
        ERInterval shiftLOB shiftHIB = shift
        shiftWidthB = shiftHIB - shiftLOB
    
    boundsAddErr errB (pLO, pHI) =
        (pLO `plusDown` (- errPoly), pHI + errPoly)
        where
        errPoly = UFB.const errB
    

instance 
    (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
    FA.ERFnApprox boxra varid ra ra (ERFnInterval fb ra)
    where
    check callerLocation f@(ERFnIntervalAny c) = f
    check callerLocation (ERFnInterval u ln c g) =
        ERFnInterval 
            (UFB.check (callerLocation ++ "upper: ") u) 
            (UFB.check (callerLocation ++ "neg lower: ") ln) 
            c g 
    domra2ranra _ = id
    ranra2domra _ = id
    setMaxDegree maxDegr (ERFnIntervalAny c) =
        ERFnIntervalAny (c { erfnMaxDegree = maxDegr } )
    setMaxDegree maxDegr (ERFnInterval u ln c g) =
        ERFnInterval 
            (UFB.reduceDegreeUp maxDegr u)
            (UFB.reduceDegreeUp maxDegr ln)
            (c { erfnMaxDegree = maxDegr } )
            g
    getRangeApprox (ERFnIntervalAny _) = RA.bottomApprox 
    getRangeApprox (ERFnInterval u ln c g) =
        UFB.raFromEndpoints u
        (
         (UFB.upperBound 10 u)
        ,
         (- (UFB.upperBound 10 ln))
        )
    scale ratio f@(ERFnIntervalAny c) = f
    scale ratio f@(ERFnInterval u ln c g) = 
        case RA.compareReals ratio 0 of
            Just GT -> 
                ERFnInterval (UFB.scaleApproxUp ratio u) (UFB.scaleApproxUp ratio ln) c g
            Just LT -> 
                ERFnInterval (UFB.scaleApproxUp (- ratio) ln) (UFB.scaleApproxUp (- ratio) u) c g
            _ -> 
                (UFA.const [ratio]) * f
    eval ptBox (ERFnIntervalAny c) = [RA.bottomApprox]
    eval ptBox (ERFnInterval u ln c g) =
        [lo RA.\/ up]
        where
        up = UFB.evalApprox ptBox u
        lo = negate $ UFB.evalApprox ptBox ln
    partialEval substitutions f@(ERFnIntervalAny c) = f
    partialEval substitutions (ERFnInterval u ln c g) =
        (ERFnInterval uP lnP c g)
        where
        uP = UFB.partialEvalApproxUp substitutions u
        lnP = UFB.partialEvalApproxUp substitutions ln

    composeThin
            f@(ERFnIntervalAny ctxt)
            substitutions =
        f
    composeThin
            f@(ERFnInterval h1 ln1 ctxt1 gl1)
            substitutions =
        (ERFnInterval h ln ctxt1 gl1)
        where
        h = UFB.composeUp maxDegree h1 ufbSubstitutions 
        ln = UFB.composeUp maxDegree ln1 ufbSubstitutions
        ufbSubstitutions = Map.map erfnUpper substitutions
        maxDegree = erfnMaxDegree ctxt1        
--        ctxt = erfnContextUnify ctxt1 ctxt2

instance 
    (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>
    UFA.ERUnitFnApprox boxra varid ra ra (ERFnInterval fb ra)
    where
    bottomApprox =
        ERFnIntervalAny erfnContextDefault
    const [val] 
        | RA.isBounded val =
            ERFnInterval
            {
                erfnUpper = fbH,
                erfnLowerNeg = fbLNeg,
                erfnContext = context,
                erfnGlobal = val
            }
        | otherwise =
            ERFnIntervalAny context 
        where
        fbH = UFB.const valH
        fbLNeg = UFB.const (negate valL)
        (valL, valH) = UFB.raEndpoints fbH val
        context = 
            erfnContextDefault
            {
                erfnCoeffGranularity = RA.getGranularity val
            }
    affine [val] coeffsSingletons
        | RA.isBounded val && (and $ map (RA.isBounded . head) $ Map.elems coeffsSingletons) =
            ERFnInterval
            {
                erfnUpper = fbH,
                erfnLowerNeg = fbLNeg,
                erfnContext = context,
                erfnGlobal = 
                    UFB.raFromEndpoints fbH
                        (valL - coeffCorr - coeffsAbsSum, 
                         valH + coeffCorr + coeffsAbsSum)
            }
        | otherwise =
            ERFnIntervalAny context
        where
        coeffs = Map.map (\[a] -> a) coeffsSingletons
        coeffGranularity =
            Map.fold max (RA.getGranularity val) (Map.map RA.getGranularity coeffs)
        coeffsMsCorrs = 
            Map.map (\(l,h) ->
                    (B.setMinGranularity coeffGranularity (l + h)/2, 
                     B.setMinGranularity coeffGranularity (h - l)/2)) $
            Map.map (UFB.raEndpoints fbH) $ coeffs
        coeffCorr = Map.fold (+) 0 $ Map.map snd coeffsMsCorrs
        coeffsAbsSum = Map.fold (+) 0 $ Map.map (abs . fst) coeffsMsCorrs
        fbH = UFB.affine (valH + coeffCorr)  (Map.map fst coeffsMsCorrs)
        fbLNeg = UFB.affine (negate (valL - coeffCorr)) (Map.map (negate . fst) coeffsMsCorrs)
        (valL, valH) = UFB.raEndpoints fbH val
        context = 
            erfnContextDefault
            {
                erfnCoeffGranularity = coeffGranularity
            }
    intersectMeasureImprovement ix vars
            f1@(ERFnIntervalAny ctxt1) 
            f2@(ERFnIntervalAny ctxt2) =
        (ERFnIntervalAny ctxt, RA.bottomApprox)
        where
        ctxt = erfnContextUnify ctxt1 ctxt2
    intersectMeasureImprovement ix vars
            f1@(ERFnIntervalAny ctxt1) 
            f2@(ERFnInterval h2 ln2 ctxt2 gl2) =
        (ERFnInterval h2 ln2 ctxt gl2, 1 / 0)
        where
        ctxt = erfnContextUnify ctxt1 ctxt2
    intersectMeasureImprovement ix vars
            f1@(ERFnInterval h1 ln1 ctxt1 gl1) 
            f2@(ERFnIntervalAny ctxt2) = 
        (ERFnInterval h1 ln1 ctxt gl1, 1)
        where
        ctxt = erfnContextUnify ctxt1 ctxt2
    intersectMeasureImprovement ix vars
            f1@(ERFnInterval h1 ln1 ctxt1 gl1) 
            f2@(ERFnInterval h2 ln2 ctxt2 gl2) =
        case RA.compareReals improvementRA 1 of
            Just LT -> (f1, 1) -- intersection made it worse, keep original
            _ ->  (intersection, improvementRA)
        where
        intersection = f1 RA./\ f2
        improvementRA 
            | 0 `RA.refines` intersectionVolume && 0 `RA.refines` f1Volume = 1
--                error $ 
--                    "ERFnInterval: intersectMeasureImprovement: inconsistent result: " 
--                    ++ show intersection
            | otherwise =
                 f1Volume / intersectionVolume
        intersectionVolume = UFA.volume vars intersection
        f1Volume = UFA.volume vars f1
        ctxt = erfnContextUnify ctxt1 ctxt2
    volume vars (ERFnIntervalAny c) = 1/0
    volume vars (ERFnInterval u ln c g) =
--        unsafePrint ("ERFnInterval: volume: result = " ++ show result) $ result
--        where
--        result =
            UFB.raFromEndpoints u $ UFB.volumeAboveZero vars (u + ln)
    integrate 
            ix (ERFnInterval u ln c g) x 
            origin (ERFnInterval uInit lnInit cInit gInit) =
--        unsafePrint
--        (
--            "ERFnInterval: integrate: " 
--            ++ "\n u = " ++ show u
--            ++ "\n ln = " ++ show ln
--            ++ "\n origin = " ++ show origin
--            ++ "\n uInit = " ++ show uInit
--            ++ "\n lnInit = " ++ show lnInit
--            ++ "\n uIuL = " ++ show uIuL
--            ++ "\n uIuU = " ++ show uIuU
--            ++ "\n uIuOriginL = " ++ show uIuOriginL
--            ++ "\n uIuOriginU = " ++ show uIuOriginU
--            ++ "\n lnIuL = " ++ show lnIuL
--            ++ "\n lnIuU = " ++ show lnIuU
--            ++ "\n lnIuOriginL = " ++ show lnIuOriginL
--            ++ "\n lnIuOriginU = " ++ show lnIuOriginU
--            ++ "\n uIov = " ++ show uIov
--            ++ "\n lnIov = " ++ show lnIov
--        )
        (ERFnInterval uIov lnIov c gIov)
        where
        -- perform raw integration of both bounds:
        (uIuL, uIuU) = 
--            mapPair (UFB.reduceDegreeDown maxDegree, UFB.reduceDegreeUp maxDegree) $ 
                UFB.integrate x u
        (lnIuL, lnIuU) = 
--            mapPair (UFB.reduceDegreeDown maxDegree, UFB.reduceDegreeUp maxDegree) $ 
                UFB.integrate x ln
        maxDegree = erfnMaxDegree c
        -- constrain the raw integrals to the origin:
        uIuOriginL = UFB.composeDown maxDegree uIuL substXOrigin
        uIuOriginU = UFB.composeUp maxDegree uIuU substXOrigin
        lnIuOriginL = UFB.composeDown maxDegree lnIuL substXOrigin
        lnIuOriginU = UFB.composeUp maxDegree lnIuU substXOrigin
        substXOrigin = Map.singleton x originUFB
        originUFB = UFB.const $ fst $ UFB.raEndpoints u origin
        -- adjust the raw integrated functions enclose the initial condition function:                        
        uIov = 
            UFB.reduceDegreeUp maxDegree $
                uIuU + uInit - uIuOriginL + (uIuOriginU - uIuOriginL)
        lnIov = 
            UFB.reduceDegreeUp maxDegree $
                lnIuU + lnInit - lnIuOriginL + (lnIuOriginU - lnIuOriginL)
        
        gIov = 
            gInit + g * ((1 - origin) RA.\/ (-1 - origin))