{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-|
    Module      :  Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate
    Description :  (testing) generating polynomials for tests
    Copyright   :  (c) 2007-2008 Michal Konecny
    License     :  BSD3

    Maintainer  :  mik@konecny.aow.cz
    Stability   :  experimental
    Portability :  portable
    
    A collection of polynomials to pick from when testing.
-}
module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Tests.Generate
where

import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Ring
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Reduce
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds
import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Enclosure

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

import Data.Number.ER.Real.DefaultRepr
import Data.Number.ER.Real.DomainBox.IntMap
import Data.Number.ER.Real.Approx.Interval
import qualified Data.Number.ER.Real.Approx as RA


import Test.QuickCheck hiding (two, three)

import qualified Data.Map as Map

{---------------------}
{----- Type synonyms for different polynomial generation distributions ----}
{---------------------}

type P = ERChebPoly (Box Int) BM

newtype PNoLimits = PNoLimits (Int, P) deriving (Show)
newtype PSize10Degree3 = PSize10Degree3 (Int, P) deriving (Show)
newtype PSize10Degree10 = PSize10Degree10 (Int, P) deriving (Show)
newtype PSize10 = PSize10 (Int, P) deriving (Show)
newtype PSize30 = PSize30 ((Int, Int), P) deriving (Show)

instance (Arbitrary PNoLimits)
    where
    arbitrary =
        elements $ map PNoLimits $ zip [0..] $ 
            polynomials1200ish id
    coarbitrary p =
        error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for polynomials"

instance (Arbitrary PSize10Degree3) 
    where
    arbitrary =
        elements $ map PSize10Degree3 $ zip [0..] $ polynomials1200ishSize10Degree3 
    coarbitrary (PSize10Degree3 p) =
        error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for polynomials"

polynomials1200ishSize10Degree3 =
    polynomials1200ish $ chplReduceTermCountUp 10 . chplReduceDegreeUp 3

instance (Arbitrary PSize10Degree10) 
    where
    arbitrary =
        elements $ map PSize10Degree10 $ zip [0..] $ 
            polynomials1200ishSize10Degree10
    coarbitrary (PSize10Degree10 p) =
        error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for polynomials"

polynomials1200ishSize10Degree10 =
    polynomials1200ish $ chplReduceTermCountUp 10 . chplReduceDegreeUp 10

instance (Arbitrary PSize10) 
    where
    arbitrary =
        elements $ map PSize10 $ zip [0..] $ polynomials1200ishSize10 
            
    coarbitrary (PSize10 p) =
        error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for polynomials"

polynomials1200ishSize10 =
    polynomials1200ish $ chplReduceTermCountUp 10
    
instance (Arbitrary PSize30) 
    where
    arbitrary =
        sized arbitrarySized
        where
        arbitrarySized n 
            | n <= 28 =
                elements $ map PSize30 $ 
                    zip (map (\n -> (0,n)) [0..]) $ 
                        polynomials200ishSize30
            | otherwise =
                elements $ map PSize30 $ 
                    zip (map (\n -> (1,n)) [0..]) $ 
                        polynomials1200ishSize30
    coarbitrary (PSize30 p) =
        error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for polynomials"

polynomials1200ishSize30 =
    polynomials1200ish $ chplReduceTermCountUp 30
    
polynomials200ishSize30 =
    polynomials200ishSmall $ chplReduceTermCountUp 30
    
data Int20 = Int20 Int deriving (Show)
    
instance (Arbitrary Int20)
    where
    arbitrary =
        do
        n <- choose (2,20)
        return $ Int20 n
    coarbitrary (Int20 n) =
        error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for EffIx20"

data Deg20Size20 = Deg20Size20 Int Int deriving (Show)
    
instance (Arbitrary Deg20Size20)
    where
    arbitrary =
        do
        maxDegree <- choose (2,20)
        maxSize <- choose (10,20)
        return $ Deg20Size20 maxDegree maxSize
    coarbitrary (Deg20Size20 maxDegree maxSize) =
        error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for Deg20Size20"

data SmallRatio = SmallRatio Int Int deriving (Show)
    
instance (Arbitrary SmallRatio)
    where
    arbitrary =
        do
        num <- choose (-1000000,1000000)
        denom <- choose (1,1000000)
        return $ SmallRatio num denom
    coarbitrary (SmallRatio num denom) =
        error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for SmallRatio"
        
        
{------------------}
{--------   Functions commonly used in tests:    ----------}
{------------------}

chplAtKeyPointsCanBeLeq ::
    (B.ERRealBase b, DomainBox box varid Int, Ord box, 
     DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb) => 
    ERChebPoly box b ->
    ERChebPoly box b ->
    Bool
chplAtKeyPointsCanBeLeq p1 p2 =
    and $ map testPoint points
    where
    points = getKeyPoints (p1 +^ p2)
    testPoint point 
        | lower1 <= upper2 =
            True
        | otherwise =
            unsafePrint
            (
                "Failure at point = " ++ (show point)
            ) $
            False
        where
        lower1 = chplEvalDown p1 point
        upper2 = chplEvalUp p2 point 
    
getKeyPoints p =
    getKeyPointsForVars $ chplGetVars p
    
getKeyPointsForVars vars =
    points
    where
    points = map DBox.fromList $ allCombinations $ map getVarPoints varDoms
    varDoms = map (\v -> (v,unitInterval)) vars
    unitInterval = ERInterval (-1) 1
    getVarPoints (var, dom) =
        (var, [domLB, domMB, domRB])
        where
        ERInterval domLB domRB = dom
        domMB = (domLB + domRB)/2

chplAtKeyPointsPointwiseBinaryDownUpConsistent ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box, 
     DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb) =>
    ((ERInterval b) -> (ERInterval b) -> (ERInterval b)) -> 
    ERChebPoly box b ->
    ERChebPoly box b ->
    (ERChebPoly box b, ERChebPoly box b) ->
    Bool
chplAtKeyPointsPointwiseBinaryDownUpConsistent raOp p1 p2 (resLow, resHigh) =
    and $ map testPoint points
    where
    points = getKeyPoints (p1 +^ p2)
    testPoint point 
        | ok = ok
        | otherwise =
            unsafePrint
            (
                "chplAtKeyPointsPointwiseBinaryDownUpConsistent failed:"
                ++ "\n point = " ++ show point
                ++ "\n raOpAtPointHigh = " ++ show raOpAtPointHigh
                ++ "\n raOpAtPointLow = " ++ show raOpAtPointLow
                ++ "\n resAtPointHigh = " ++ show resAtPointHigh
                ++ "\n resAtPointLow = " ++ show resAtPointLow
            )
            ok
        where
        ok = 
            raOpAtPointLow <= resAtPointHigh
            &&
            raOpAtPointHigh >= resAtPointLow
        resAtPointLow = chplEvalDown resLow point
        resAtPointHigh = chplEvalUp resHigh point
        raOpAtPoint@(ERInterval raOpAtPointLow raOpAtPointHigh) = 
            raOp p1AtPoint p2AtPoint 
        p1AtPoint = ERInterval p1AtPointLow p1AtPointHigh
        (p1AtPointLow, p1AtPointHigh) = chplEval p1 point
        p2AtPoint = ERInterval p2AtPointLow p2AtPointHigh
        (p2AtPointLow, p2AtPointHigh) = chplEval p2 point

makeThickEncl maxDegree maxSize p1 p2 =
    (pMax q1Neg q2Neg, pMax q1 q2)
    where
    q1Neg = chplNeg q1
    q2Neg = chplNeg q2
    q1 = p1 +^ p2Mp1ScaledDown
    q2 = p1 -^ p2Mp1ScaledDown
    p2Mp1ScaledDown =
        chplScaleUp (10/sizeB) p2Mp1
        where
        sizeB = max (abs upperB) (abs lowerB)
        (lowerB, upperB) = chplBounds 10 p2Mp1
        p2Mp1 = p2 -^ p1
    pMax = chplMaxUp maxDegree maxSize
    
makeParalEncl p num denom =
--    unsafePrintReturn
--    (
--        "makeThinEncl: result = "
--    )
    (pNeg, p +^ cP)
    where
    pNeg = chplNeg p
    cP = chplConst cB
    cB = abs $ numB / (1000 * denomB)
    numB = fromInteger $ toInteger num
    denomB = fromInteger $ toInteger denom
    
enclRestrictRange ix (Nothing, Nothing) pEncl = (True, pEncl)
enclRestrictRange ix (maybeLower, maybeUpper) preEncl =
    (succeeded, pEncl)
    where
    succeeded = lowerSucceeded && upperSucceeded
    lowerSucceeded =
        case maybeLower of
            Nothing -> True
            Just lower -> pLowerBound > lower 
    upperSucceeded =
        case maybeUpper of
            Nothing -> True
            Just upper -> pUpperBound < upper
    (pLowerBound, pUpperBound) = enclBounds ix pEncl
    pEncl =
        case (maybeLower, maybeUpper) of
            (Just lowerB, Nothing) ->
                case lowerB <= preLowerBoundB of
                    True -> preEncl -- enclosure already in the range
                    False -> -- a shift needed to get above the lower bound
                        enclAddConst (lowerB - preLowerBoundB + sepB) preEncl
            (Nothing, Just upperB) ->
                case preUpperBoundB <= upperB of
                    True -> preEncl -- enclosure already in the range
                    False -> -- a shift needed to get below the upper bound
                        enclAddConst (upperB - preUpperBoundB - sepB) preEncl
            (Just lowerB, Just upperB) ->
                case lowerB <= preLowerBoundB && preUpperBoundB <= upperB of
                    True -> preEncl -- enclosure already in the range
                    _ -> 
                        case preWidthB + sepB <= widthB of
                            True -> -- no scaling needed, only shifting by a constant to the centre of the range
                                enclAddConst 
                                    (lowerB - preLowerBoundB + (preWidthB - widthB)/2) 
                                    preEncl
                            _ -> -- full affine transformation needed
                                enclAddConst
                                    (lowerB + sepB) $
                                    enclScaleNonneg -- scale preEncl so that it fits inside the range
                                        (widthB / saferPreWidthB) $
                                        enclAddConst -- shift preEncl so that it is non-negative and as close to 0 as safely possible
                                            (sepB - preLowerBoundB)
                                            preEncl
                where 
                widthB = upperB - lowerB
                saferPreWidthB = preWidthB + 2 * sepB
    sepB = preWidthB / 1000000
    preWidthB = preUpperBoundB - preLowerBoundB
    (preLowerBoundB, preUpperBoundB) = enclBounds ix preEncl
    
    

enclAtKeyPointsPointwiseBinaryDownUpConsistent ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box, 
     DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb, Show testId) =>
    String {-^ report file name -} ->
    testId {-^ item to identify the random input given to the test -} ->
    ((ERInterval b) -> (ERInterval b) -> (ERInterval b))
        {-^ this real approx operation has to return an inner approximation of the exact result set, 
            ie each number that the approximation supports is in the maximal extension -} ->
    (ERChebPoly box b, ERChebPoly box b) {-^ enclosure of argument 1 -} ->
    (ERChebPoly box b, ERChebPoly box b) {-^ enclosure of argument 2 -} ->
    (ERChebPoly box b, ERChebPoly box b) {-^ alleged enclosure of result -} ->
    Bool
enclAtKeyPointsPointwiseBinaryDownUpConsistent
        reportFileName testId
        raOpInner 
        p1Encl@(p1LowNeg, p1High) p2Encl@(p2LowNeg, p2High) resEncl =
    and $ map testPoint points
    where
    points = getKeyPoints (p1High +^ p2High +^ p1LowNeg +^ p2LowNeg)
    testPoint point 
        | result =
            unsafeReport reportFileName
            (
                show $ 
                    (testId, point, p1OpInnerP2AtPoint, resAtPoint)
            ) 
            result
        | otherwise = 
            unsafePrint
            (
                "enclAtKeyPointsPointwiseBinaryDownUpConsistent failed"
                ++ "\n point = " ++ show point
                ++ "\n p1AtPoint = " ++ show p1AtPoint
                ++ "\n p2AtPoint = " ++ show p2AtPoint
                ++ "\n p1OpInnerP2AtPoint = " ++ show p1OpInnerP2AtPoint
                ++ "\n resAtPoint = " ++ show resAtPoint
            ) $
            result
        where
        result = p1OpInnerP2AtPoint `RA.refines` resAtPoint
        p1OpInnerP2AtPoint = p1AtPoint `raOpInner` p2AtPoint
        resAtPoint = enclEval resEncl point
--        resAtPoint = p1OpInnerP2AtPoint -- for dummy testing that never <<loop>>s
        p1AtPoint = normaliseERInterval $ enclEvalInner p1Encl point
        p2AtPoint = normaliseERInterval $ enclEvalInner p2Encl point

enclAtKeyPointsPointwiseUnaryDownUpConsistent ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box, 
     DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb, Show testId) =>
    String {-^ report file name -} ->
    testId {-^ item to identify the random input given to the test -} ->
    ((ERInterval b) -> (ERInterval b))
        {-^ this real approx operation has to return an inner approximation of the exact result set, 
            ie each number that the approximation supports is in the maximal extension -} ->
    (ERChebPoly box b, ERChebPoly box b) {-^ enclosure of argument -} ->
    (ERChebPoly box b, ERChebPoly box b) {-^ alleged enclosure of result -} ->
    Bool
enclAtKeyPointsPointwiseUnaryDownUpConsistent
        reportFileName testId
        raOpInner 
        pEncl@(pLowNeg, pHigh) resEncl =
    and $ map testPoint points
    where
    points = getKeyPoints (pHigh +^ pLowNeg)
    testPoint point 
        | result =
            unsafeReport reportFileName
            (
                show $ 
                    (testId, point, opInnerPAtPoint, resAtPoint)
            )
            result 
        | otherwise = 
            unsafePrint
            (
                "enclAtKeyPointsPointwiseUnaryDownUpConsistent failed"
                ++ "\n point = " ++ show point
                ++ "\n pAtPoint = " ++ show pAtPoint
                ++ "\n opInnerPAtPoint = " ++ show opInnerPAtPoint
                ++ "\n resAtPoint = " ++ show resAtPoint
            ) $
            result
        where
        result = opInnerPAtPoint `RA.refines` resAtPoint
        opInnerPAtPoint = raOpInner pAtPoint
        resAtPoint = enclEval resEncl point
        pAtPoint = 
--            normaliseERInterval $ 
            enclEvalInner pEncl point


enclAtKeyPointsConsistent ::
    (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box, 
     DomainBoxMappable boxb boxbb varid b [ERInterval b], Show boxb, Show testId) =>
    String {-^ report file name -} ->
    testId {-^ item to identify the random input given to the test -} ->
    (boxb -> (ERInterval b))
        {-^ this operation has to return an inner approximation of the exact result set, 
            ie each number that the approximation supports is a solution in the maximal extension -} ->
    [varid] {-^ variables to test over -} ->
    (ERChebPoly box b, ERChebPoly box b) {-^ alleged enclosure of result -} ->
    Bool
enclAtKeyPointsConsistent
        reportFileName testId
        opInner allVars resEncl@(resLowNeg, resHigh) =
    and $ map testPoint points
    where
    points = getKeyPointsForVars allVars
    testPoint point 
        | result =
            unsafeReport reportFileName
            (
                show $ 
                    (testId, point, opInnerAtPoint, resAtPoint)
            )
            result 
        | otherwise = 
            unsafePrint
            (
                "enclAtKeyPointsConsistent failed"
                ++ "\n point = " ++ show point
                ++ "\n opInnerAtPoint = " ++ show opInnerAtPoint
                ++ "\n resAtPoint = " ++ show resAtPoint
            ) $
            result
        where
        result = opInnerAtPoint `RA.refines` resAtPoint
        opInnerAtPoint = opInner point
        resAtPoint = enclEval resEncl point


{------------------}
{--------   A diverse collection of polynomials to pick from:    ----------}
{------------------}

type E = (P,P)

vars :: [P]
vars = map chplVar [0..7]

varsE :: [E]
varsE = map (\p -> (chplNeg p, p)) vars

x0 = vars !! 0
x1 = vars !! 1
x2 = vars !! 2
x3 = vars !! 3
x4 = vars !! 4

x0E = varsE !! 0
x1E = varsE !! 1
x2E = varsE !! 2
x3E = varsE !! 3
x4E = varsE !! 4

one :: P
[mone, one, two, three, seven, thousand, million, tiny, huge] = 
    map chplConst 
    [-1,1,2,3,7,1000,1000000,10^^(-200),10^^200]

oneE :: E
[moneE, oneE, twoE, threeE, sevenE, thousandE, millionE, tinyE, hugeE] = 
    map (\ c -> (chplConst (-c), chplConst c))
    [-1,1,2,3,7,1000,1000000,10^^(-200),10^^200]

polynomials1200ish rdc =
    concat $ map (powers10 rdc) $
    concat $ map addConsts3 $
    concat $ map multConsts3 $
    polyBase13
    
polynomials200ish rdc =
    concat $ map (powers4 rdc) $
    concat $ map addConsts3 $
    concat $ map multConsts3 $
    polyBase5
    
polynomials40ish rdc =
    concat $ map (powers2 rdc) $
    concat $ map addConsts2 $
    concat $ map multConsts2 $
    polyBase5
    
polynomials200ishSmall rdc =
    concat $ map (powers4Small rdc) $
    concat $ map addConsts3 $
    concat $ map multConsts3 $
    polyBase5
    
polynomials40ishSmall rdc =
    concat $ map (powers2Small rdc) $
    concat $ map addConsts2 $
    concat $ map multConsts2 $
    polyBase5
    

polyBase5 =
        [
         (two *^ x0) +^ x1
        ,
         (seven *^ x0) -^ x1
        ,
         (tiny *^ x0) +^ x1
        ,
         x0 -^ x1 *^ x2
        ,
         x0 -^ x1 +^ x2 -^ x3 +^ x4
        ]
    
polyBase13 =
        [
         x0
        ,
         x0 +^ x1
        ,
         x0 -^ x1
        ,
         (two *^ x0) +^ x1
        ,
         (two *^ x0) -^ x1
        ,
         (seven *^ x0) +^ x1
        ,
         (seven *^ x0) -^ x1
        ,
         (tiny *^ x0) +^ x1
        ,
         (tiny *^ x0) -^ x1
        ,
         x0 -^ x1 +^ x2
        ,
         x0 -^ x1 *^ x2
        ,
         x0 +^ x1 +^ x2 +^ x3 +^ x4
        ,
         x0 -^ x1 +^ x2 -^ x3 +^ x4
        ]
    
powersAll rdc p =
    powersAux [p, rdc $ p *^ p]
    where
    powersAux (pNHalfM1 : pNHalf : rest) = 
        pNHalfM1 : (powersAux $ (pNHalf : rest) ++ [pNM1, pN])
        where
        pNM1 = rdc $ pNHalf *^ pNHalfM1
        pN = rdc $ pNHalf *^ pNHalf

powersForExps rdc p exponents =
    map pw exponents
    where
    pw n = pws !! (n - 1)
    pws = powersAll rdc p

powers10 rdc p =
    powersForExps rdc p [1..10]

powers4 rdc p =
    powersForExps rdc p [1,3,5,7]
    
powers4Small rdc p =
    powersForExps rdc p [1,2,3,5]
    
powers2 rdc p =
    powersForExps rdc p [1,7]
    
powers2Small rdc p =
    powersForExps rdc p [1,3]
    
addConsts3 p =
    [p +^ one, p +^ three, p +^ seven]

multConsts3 p =
    [p *^ two, p *^ three, p *^ seven]
    
addConsts2 p =
    [p +^ one, p +^ three]

multConsts2 p =
    [p *^ two, p *^ seven]