{-# LANGUAGE CPP #-}
{-# OPTIONS_GHC -fno-warn-missing-methods #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE UndecidableInstances   #-}
{-# LANGUAGE FlexibleInstances   #-}
{-# LANGUAGE FlexibleContexts   #-}
{-# LANGUAGE DeriveDataTypeable   #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-|
    Module      :  Data.Number.ER.RnToRm.Approx.DomTransl
    Description :  enclosures translated from [-1,1]^n to another domain 
    Copyright   :  (c) Michal Konecny
    License     :  BSD3

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

    Datatype translating enclosures from @[-1,1]^n@ to any compact
    interval in @R^n@ with non-empty interior.
-}
-- #define ASSUME_DOMAINS_COMPATIBLE
module Data.Number.ER.RnToRm.Approx.DomTransl 
(
    ERFnDomTranslApprox(..), DomTransl(..)
)
where

import qualified Data.Number.ER.RnToRm.Approx as FA
import qualified Data.Number.ER.RnToRm.UnitDom.Approx as UFA
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(..), DomainIntBox, DomainBoxMappable)
import Data.Number.ER.BasicTypes
import Data.Number.ER.Misc

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

import qualified Data.Map as Map

{-|
    Datatype translating enclosures from @[-1,1]^n@ to any compact
    interval in @R^n@ with non-empty interior
    using a bunch of linear maps, one for each dimension. 
-}
data ERFnDomTranslApprox dtrbox varid ufa ira =
    ERFnDomTranslApprox
    {
        erfnUnitApprox :: ufa,
        erfnDomTransl :: dtrbox
    }
    deriving (Typeable, Data)

instance (Binary a, Binary b, Binary c, Binary d) => Binary (ERFnDomTranslApprox a b c d) where
  put (ERFnDomTranslApprox a b) = put a >> put b
  get = get >>= \a -> get >>= \b -> return (ERFnDomTranslApprox a b)
    
{-| 
    The canonical translation of 
    any compact non-empty and non-singleton interval in @R@
    to and from the unit interval @[-1,1]@.
    
    This structure holds the two coefficients for both
    linear mappings.
-}
data DomTransl ira =
    DomTransl
    {
        dtrDom :: ira {-^ the interval being mapped -},
        dtrFromUnitSlope :: ira,
        dtrFromUnitConst :: ira,
        dtrToUnitSlope :: ira,
        dtrToUnitConst :: ira
    }
    deriving (Typeable, Data)
    
instance (Binary a) => Binary (DomTransl a) where
  put (DomTransl a b c d e) = put a >> put b >> put c >> put d >> put e
  get = get >>= \a -> get >>= \b -> get >>= \c -> get >>= \d -> get >>= \e -> return (DomTransl a b c d e)
    
instance
    (RA.ERIntApprox domra) =>
    Eq (DomTransl domra)
    where
    (DomTransl _ _ _ _ dom1) == (DomTransl _ _ _ _ dom2) =
        RA.equalApprox dom1 dom2
    
instance
    (RA.ERIntApprox domra) =>
    Show (DomTransl domra)
    where
    show (DomTransl fromA fromB toA toB dom) =
        "DomTransl\n" ++   
        "  dom = " ++ show dom ++ "\n" ++
        "  fromUnit = " ++ show fromA ++ " * x + " ++ show fromB ++ "\n" ++
        "  toUnit = " ++ show toA ++ " * x + " ++ show toB ++ "\n"

dtrIdentity ::
    (RA.ERIntApprox ira) =>
    DomTransl ira
dtrIdentity =
    makeDomTransl ((-1) RA.\/ 1)
    
dtrBToDomB dtrB =
    DBox.map dtrDom dtrB
    
makeDomTransl ::
    (RA.ERIntApprox ira) =>
    ira ->
    DomTransl ira
makeDomTransl dom 
    | domSuitable =
        DomTransl
        {
            dtrFromUnitSlope = dHMdL / 2,
            dtrFromUnitConst = dHPdL / 2,
            dtrToUnitSlope = 2 / dHMdLgr,
            dtrToUnitConst = - dHPdL / dHMdLgr,
            dtrDom = dom
        }
    | otherwise =
        error $ 
            "DomTranslApprox: makeDomTransl: cannot make a translation to domain " 
            ++ show dom
    where
    domSuitable = RA.isBounded dom && (not $ RA.isExact dom)
    (dL, dH) = RA.bounds dom
    dHPdL = dH + dL
    dHMdL = dH - dL
    dHMdLgr = RA.setMinGranularity 100 dHMdL
--    fromUnit x = (x * (dHMdL) + dHPdL) / 2 
--    toUnit y = (2 * y - dHPdL) / dHMdL

dtrToUnit domTransl x = a * x + b
    where
    a = dtrToUnitSlope domTransl
    b = dtrToUnitConst domTransl
    
dtrFromUnit domTransl x = a * x + b
    where
    a = dtrFromUnitSlope domTransl
    b = dtrFromUnitConst domTransl

domToUnit ::
--    (DomainIntBox dbox varid ira, 
--     DomainIntBox dtrbox varid (DomTransl ira)) =>
    (DomainBoxMappable dbox dtrbox varid ira (DomTransl ira),
     Num ira) => 
    dtrbox -> dbox -> dbox
domToUnit dtrB domBox =
    DBox.intersectionWith (\d dtr -> dtrToUnit dtr d) domBox dtrB

#ifdef ASSUME_DOMAINS_COMPATIBLE

dtrsCompatible _ _ = True

dtrUnion msg dtr1 dtr2 = dtr1

#else    
dtrsCompatible dtr1 dtr2 =
    foldl (&&) True $ map snd $ 
        DBox.zipWith eqDomains dtr1 dtr2
    where
    eqDomains d1 d2 =
        d1L == d2L && d1U == d2U
        where
        (d1L, d1U) = RA.bounds $ dtrDom d1
        (d2L, d2U) = RA.bounds $ dtrDom d2

dtrUnion msg dtr1 dtr2 
    | dtrsCompatible dtr1 dtr2 = 
        DBox.union dtr1 dtr2
    | otherwise = error msg

#endif

dtrBShow dtrs =
    concatWith "," $
    map showOne $ DBox.toList dtrs
    where
    showOne (var, dtr) =
        showVar var ++ " in " ++ show (dtrDom dtr)
        
    
instance 
    (UFA.ERUnitFnApprox box varid domra ranra ufa, 
     DomainBoxMappable dtrbox box varid (DomTransl domra) domra) =>
    Show (ERFnDomTranslApprox dtrbox varid ufa domra)
    where
    show (ERFnDomTranslApprox ufa dtrB) =
        "\nERFnDomTranslApprox" ++
        show ufaDom ++
--        show ufa ++
        "\n dom = [" ++
        (concatWith ", " $ map showVarDom $ DBox.toList $ dtrBToDomB dtrB)
        ++ "]"
        where
        ufaDom =
            FA.composeThin ufa $ 
                Map.fromAscList $ 
                    map mkToUnitUFA $ 
                        DBox.toAscList dtrB
--        gr = 20 + (RA.getGranularity ufa)
        mkToUnitUFA (var, tr) =
            (var, UFA.affine [co] (Map.singleton var [sl]))
            where
            sl = FA.domra2ranra ufa $ dtrToUnitSlope tr
            co = FA.domra2ranra ufa $ dtrToUnitConst tr
        showVarDom (varID, varDom) =
            showVar varID ++ " -> " ++ show varDom 
instance 
    (UFA.ERUnitFnApprox box varid domra ranra ufa, 
     Eq dtrbox) =>
    Eq (ERFnDomTranslApprox dtrbox varid ufa domra)
    where
    (ERFnDomTranslApprox ufa1 dtrB1) == (ERFnDomTranslApprox ufa2 dtrB2) =
        ufa1 == ufa2 && dtrB1 == dtrB2

instance 
    (UFA.ERUnitFnApprox box varid domra ranra ufa, Ord ufa
    , Eq dtrbox) =>
    Ord (ERFnDomTranslApprox dtrbox varid ufa domra)
    where
    compare (ERFnDomTranslApprox ufa1 dtrB1) (ERFnDomTranslApprox ufa2 dtrB2)
        | dtrB1 == dtrB2 =
                compare ufa1 ufa2
        | otherwise =
            error "DomTransl: compare: incompatible domains"  

instance
    (UFA.ERUnitFnApprox box varid domra ranra ufa, 
     DomainBoxMappable dtrbox box varid (DomTransl domra) domra, Eq dtrbox) =>
    Num (ERFnDomTranslApprox dtrbox varid ufa domra)
    where
    fromInteger n = ERFnDomTranslApprox (fromInteger n) DBox.noinfo
    negate (ERFnDomTranslApprox ufa dtrB) =
        (ERFnDomTranslApprox (negate ufa) dtrB)
    (ERFnDomTranslApprox ufa1 dtr1) + (ERFnDomTranslApprox ufa2 dtr2) =
        ERFnDomTranslApprox (ufa1 + ufa2) (dtrUnion msg dtr1 dtr2)
        where
        msg = "DomTransl: cannot add approximations with incompatible domains"
    (ERFnDomTranslApprox ufa1 dtr1) * (ERFnDomTranslApprox ufa2 dtr2) = 
        ERFnDomTranslApprox (ufa1 * ufa2) (dtrUnion msg dtr1 dtr2)
        where
        msg = "DomTransl: cannot multiply approximations with incompatible domains"
        
instance 
    (UFA.ERUnitFnApprox box varid domra ranra ufa
    , DomainBoxMappable dtrbox box varid (DomTransl domra) domra, Eq dtrbox) =>
    Fractional (ERFnDomTranslApprox dtrbox varid ufa domra)
    where
    fromRational r = ERFnDomTranslApprox (fromRational r) DBox.noinfo
    recip (ERFnDomTranslApprox ufa dtrB) =
        ERFnDomTranslApprox (recip ufa) dtrB

instance 
    (UFA.ERUnitFnApprox box varid domra ranra ufa
    , DomainBoxMappable dtrbox box varid (DomTransl domra) domra, Eq dtrbox) =>
    RA.ERApprox (ERFnDomTranslApprox dtrbox varid ufa domra)
    where
    initialiseBaseArithmetic _ =
    	RA.initialiseBaseArithmetic (0 :: ufa)
    getGranularity (ERFnDomTranslApprox ufa dtrB) =
        RA.getGranularity ufa
    setGranularity gran (ERFnDomTranslApprox ufa dtrB) =
        ERFnDomTranslApprox (RA.setGranularity gran ufa) dtrB
    setMinGranularity gran (ERFnDomTranslApprox ufa dtrB) =
        ERFnDomTranslApprox (RA.setMinGranularity gran ufa) dtrB
    (ERFnDomTranslApprox ufa1 dtrB1) /\ (ERFnDomTranslApprox ufa2 dtrB2) =
        ERFnDomTranslApprox (ufa1 RA./\ ufa2) (dtrUnion msg dtrB1 dtrB2)
        where
        msg = "DomTransl: cannot intersect approximations with incompatible domains"
    intersectMeasureImprovement ix 
            (ERFnDomTranslApprox ufa1 dtrB1) 
            (ERFnDomTranslApprox ufa2 dtrB2) = 
        (ERFnDomTranslApprox ufaIsect dtrB,
         ERFnDomTranslApprox ufaImpr dtrB)
        where
        (ufaIsect, raImpr) = UFA.intersectMeasureImprovement ix vars ufa1 ufa2
        ufaImpr = UFA.const [raImpr]
        dtrB = dtrUnion msg dtrB1 dtrB2
        msg = "DomTransl: cannot intersect approximations with incompatible domains"
        vars = DBox.keys dtrB
    leqReals fa1 fa2 =
        RA.leqReals (erfnUnitApprox fa1) (erfnUnitApprox fa2)


instance
    (UFA.ERUnitFnApprox box varid domra ranra ufa, RA.ERIntApprox ufa 
    , DomainBoxMappable dtrbox box varid (DomTransl domra) domra, Eq dtrbox) =>
    RA.ERIntApprox (ERFnDomTranslApprox dtrbox varid ufa domra)
    where
--    doubleBounds = :: ira -> (Double, Double) 
--    floatBounds :: ira -> (Float, Float)
--    integerBounds :: ira -> (ExtendedInteger, ExtendedInteger)
    bisectDomain maybePt (ERFnDomTranslApprox ufa dtrB) =
        (ERFnDomTranslApprox ufa1 dtrB,
         ERFnDomTranslApprox ufa2 dtrB)
         where
         (ufa1, ufa2) = RA.bisectDomain (fmap erfnUnitApprox maybePt) ufa
    bounds (ERFnDomTranslApprox ufa dtrB) =
        (ERFnDomTranslApprox ufa1 dtrB,
         ERFnDomTranslApprox ufa2 dtrB)
         where
         (ufa1, ufa2) = RA.bounds ufa
    (ERFnDomTranslApprox ufa1 dtrB1) \/ (ERFnDomTranslApprox ufa2 dtrB2) =
        ERFnDomTranslApprox (ufa1 RA.\/ ufa2) (dtrUnion msg dtrB1 dtrB2)
        where
        msg = "DomTransl: cannot intersect approximations with incompatible domains"

instance
    (UFA.ERUnitFnApprox box varid domra ranra ufa, RAEL.ERApproxElementary ufa
    , DomainBoxMappable dtrbox box varid (DomTransl domra) domra, Eq dtrbox) =>
    RAEL.ERApproxElementary (ERFnDomTranslApprox dtrbox varid ufa domra)
    where
    abs ix (ERFnDomTranslApprox ufa dtrB) =
        ERFnDomTranslApprox (RAEL.abs ix ufa) dtrB
    exp ix (ERFnDomTranslApprox ufa dtrB) =
        ERFnDomTranslApprox (RAEL.exp ix ufa) dtrB
    log ix (ERFnDomTranslApprox ufa dtrB) =
        ERFnDomTranslApprox (RAEL.log ix ufa) dtrB
    sin ix (ERFnDomTranslApprox ufa dtrB) =
        ERFnDomTranslApprox (RAEL.sin ix ufa) dtrB
    cos ix (ERFnDomTranslApprox ufa dtrB) =
        ERFnDomTranslApprox (RAEL.cos ix ufa) dtrB
    atan ix (ERFnDomTranslApprox ufa dtrB) =
        ERFnDomTranslApprox (RAEL.atan ix ufa) dtrB

instance 
    (UFA.ERUnitFnApprox box varid domra ranra ufa, 
     DomainIntBox box varid domra, 
     DomainBoxMappable dtrbox box varid (DomTransl domra) domra, 
     DomainBoxMappable box dtrbox varid domra (DomTransl domra), 
     Eq dtrbox) =>
    FA.ERFnApprox box varid domra ranra (ERFnDomTranslApprox dtrbox varid ufa domra)
    where
    check prgLocation (ERFnDomTranslApprox ufa dtrB) =
        ERFnDomTranslApprox (FA.check (prgLocation ++ dtrBShow dtrB ++ ": ") ufa) dtrB
    domra2ranra fa d =
        FA.domra2ranra (erfnUnitApprox fa) d
    ranra2domra fa r =
        FA.ranra2domra (erfnUnitApprox fa) r
    setMaxDegree maxDegree (ERFnDomTranslApprox ufa dtrB) =
        ERFnDomTranslApprox (FA.setMaxDegree maxDegree ufa) dtrB
    volume (ERFnDomTranslApprox ufa dtrB) =
        DBox.fold 
            (\tr vol -> vol * (FA.domra2ranra ufa $ dtrFromUnitSlope tr)) 
            (UFA.volume vars ufa) dtrB 
        where
        vars = DBox.keys dtrB
    scale ratio (ERFnDomTranslApprox ufa dtrB) =
        (ERFnDomTranslApprox (FA.scale ratio ufa) dtrB)
    partialIntersect ix substitutions f1 f2 
        | insideSubstitutions = f1 RA./\ f2
        | otherwise = f1
        where
        insideSubstitutions =
            and $ map snd $
                DBox.zipWith (RA.refines) dom1 substitutions
        dom1 = FA.dom f1
    eval ptBox (ERFnDomTranslApprox ufa dtrB) =
        FA.eval (domToUnit dtrB ptBox) ufa
    partialEval substitutions (ERFnDomTranslApprox ufa dtrB) =
        (ERFnDomTranslApprox (FA.partialEval (domToUnit dtrB substitutions) ufa) dtrBNoVars)
        where
        dtrBNoVars =
            DBox.difference dtrB substitutions

--instance 
--    (UFA.ERUnitFnApprox box varid domra ranra ufa, 
--     DomainIntBox box varid domra, 
--     VariableID varid) =>
--    UFA.ERUnitFnApprox box varid domra ranra (ERFnDomTranslApprox dtrbox varid ufa domra)
--    where
--    const vals =
--        ERFnDomTranslApprox
--        {
--            erfnUnitApprox = UFA.const vals,
--            erfnDomTransl = Map.empty
--        }
--    affine c coeffs =
--        ERFnDomTranslApprox
--        {
--            erfnUnitApprox = UFA.affine c coeffs,
--            erfnDomTransl = Map.map (const dtrIdentity) coeffs
--        }

instance 
    (UFA.ERUnitFnApprox box varid domra ranra ufa,
     DomainIntBox box varid domra,
     DomainBoxMappable dtrbox box varid (DomTransl domra) domra, 
     DomainBoxMappable box dtrbox varid domra (DomTransl domra), 
     Eq dtrbox) =>
    FA.ERFnDomApprox box varid domra ranra (ERFnDomTranslApprox dtrbox varid ufa domra)
    where
    dom (ERFnDomTranslApprox ufa dtrB) = dtrBToDomB dtrB
    bottomApprox domB tupleSize 
        | tupleSize == 1 =
            ERFnDomTranslApprox
            {
                erfnUnitApprox = UFA.bottomApprox,
                erfnDomTransl = DBox.map makeDomTransl domB
            }
    const domB vals =
        ERFnDomTranslApprox
        {
            erfnUnitApprox = UFA.const vals,
            erfnDomTransl = DBox.map makeDomTransl domB
        }
    proj domB i =
        ERFnDomTranslApprox
        {
            erfnUnitApprox = ufa,
            erfnDomTransl = domTransls
        }
        where
        domTransls = DBox.map makeDomTransl domB
        idomTransl = DBox.lookup "ERFnDomTranslApprox: ERFnDomApprox: proj: " i domTransls
        sl = FA.domra2ranra ufa $ dtrFromUnitSlope idomTransl
        co = FA.domra2ranra ufa $ dtrFromUnitConst idomTransl
        ufa = UFA.affine [co] (Map.singleton i [sl])
    -- split the function by its domain into two halves:
    bisect var maybePt f@(ERFnDomTranslApprox ufa dtrB)
        | varAbsent =
            (f, f)
        | ptOK = 
            (ERFnDomTranslApprox ufaLeft dtrLeft, 
             ERFnDomTranslApprox ufaRight dtrRight)
        | otherwise =
            error $
                "DomTransl: faBisect: bisection point " ++ show pt ++
                " is not exact " ++
                "(var = " ++ showVar var ++ ")" ++ 
                "(domain = " ++ show dom ++ ")"
        where
        (pt, ptOK) = 
            case maybePt of
                Just pt -> (pt, RA.isExact pt)
                Nothing -> (domM, True) 
        (domL, domM, domR, domGran) = RA.exactMiddle dom
        varAbsent = DBox.notMember var dtrB
        dom = 
            dtrDom $ DBox.lookup errMsg var dtrB
            where
            errMsg =
                "ERFnDomTranslApprox: FA.bisect: var " ++ showVar var 
                ++ " not in the domain of " ++ show f
        ufaLeft = FA.composeThin ufa $ Map.singleton var toLeft 
        ufaRight = FA.composeThin ufa $ Map.singleton var toRight
        dtrLeft = DBox.insert var (makeDomTransl domLeft) dtrB 
        dtrRight = DBox.insert var (makeDomTransl domRight) dtrB
        domLeft = domL RA.\/ pt
        domRight = pt RA.\/ domR
        toLeft =
            UFA.affine [midLeft] (Map.singleton var [slopeLeft])
        toRight =
            UFA.affine [midRight] (Map.singleton var [slopeRight])
        (midLeft, slopeLeft, midRight, slopeRight) =
            getExactTransforms initGran
        initGran =
            max domGran $ RA.getGranularity pt
        getExactTransforms gran 
            | and $ map RA.isExact [midLeft, slopeLeft, midRight, slopeRight] =
                 (midLeft, slopeLeft, midRight, slopeRight)
            | otherwise = getExactTransforms (gran + 1)
            where
            midLeft = slopeLeft - 1
            midRight = 1 - slopeRight
            slopeLeft = sizeLeft / size        
            slopeRight = sizeRight / size        
            size = domRgr - domLgr
            sizeLeft = ptGr - domLgr
            sizeRight = domRgr - ptGr
            domRgr = RA.setMinGranularity gran $ FA.domra2ranra ufa domR
            domLgr = RA.setMinGranularity gran $ FA.domra2ranra ufa domL
            ptGr = RA.setMinGranularity gran $ FA.domra2ranra ufa pt
    integrate
            ix fD@(ERFnDomTranslApprox ufaD dtrBD) x integdomBox
            origin (ERFnDomTranslApprox ufaInit dtrBInit) =
        ERFnDomTranslApprox ufaI dtrBD
        where
        ufaI =
            UFA.integrate
                ix ufaDadj x 
                (dtrToUnit trX origin) 
                ufaInit
        ufaDadj = 
            FA.scale (FA.domra2ranra ufaD $ dtrFromUnitSlope trX) $
            ufaD
        trX = 
            DBox.findWithDefault err x dtrBD
        err = 
            error $
                "DomTransl: faIntegrate: variable " ++ showVar x ++
                " not in the domain of the function " ++ show fD