{-# 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 <>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]