aern2-mfun (empty) → 0.2.9.0
raw patch · 14 files changed
+2118/−0 lines, 14 filesdep +aern2-fundep +aern2-mfundep +aern2-mp
Dependencies added: aern2-fun, aern2-mfun, aern2-mp, base, collect-errors, mixed-types-num, vector
Files
- aern2-mfun.cabal +91/−0
- bench/OptimisationBenchmark.hs +30/−0
- src/AERN2/AD/Differential.hs +11/−0
- src/AERN2/AD/GenericOperations.hs +220/−0
- src/AERN2/AD/MPBallOperations.hs +223/−0
- src/AERN2/AD/Type.hs +46/−0
- src/AERN2/BoxFun/Box.hs +150/−0
- src/AERN2/BoxFun/Optimisation.hs +319/−0
- src/AERN2/BoxFun/TestFunctions.hs +425/−0
- src/AERN2/BoxFun/Type.hs +155/−0
- src/AERN2/Linear/Matrix/Inverse.hs +126/−0
- src/AERN2/Linear/Matrix/Type.hs +132/−0
- src/AERN2/Linear/Vector/Type.hs +176/−0
- src/AERN2/Util/Util.hs +14/−0
+ aern2-mfun.cabal view
@@ -0,0 +1,91 @@+cabal-version: 1.12++-- This file has been generated from package.yaml by hpack version 0.34.4.+--+-- see: https://github.com/sol/hpack++name: aern2-mfun+version: 0.2.9.0+synopsis: Multi-variate real function optimisation and proving+description: Please see the README on GitHub at <https://github.com/michalkonecny/aern2/#readme>+category: Math+homepage: https://github.com/michalkonecny/aern2#readme+bug-reports: https://github.com/michalkonecny/aern2/issues+author: Eike Neumann, Junaid Rasheed, Michal Konecny+maintainer: mikkonecny@gmail.com+copyright: 2019-2022 Eike Neumann, Junaid Rasheed, Michal Konecny+license: BSD3+build-type: Simple++source-repository head+ type: git+ location: https://github.com/michalkonecny/aern2++library+ exposed-modules:+ AERN2.AD.Differential+ AERN2.AD.GenericOperations+ AERN2.AD.MPBallOperations+ AERN2.AD.Type+ AERN2.BoxFun.Box+ AERN2.BoxFun.Optimisation+ AERN2.BoxFun.TestFunctions+ AERN2.BoxFun.Type+ AERN2.Linear.Matrix.Inverse+ AERN2.Linear.Matrix.Type+ AERN2.Linear.Vector.Type+ AERN2.Util.Util+ other-modules:+ Paths_aern2_mfun+ hs-source-dirs:+ src+ default-extensions:+ RebindableSyntax,+ ScopedTypeVariables,+ TypeFamilies,+ TypeOperators,+ ConstraintKinds,+ DefaultSignatures,+ MultiParamTypeClasses,+ FlexibleContexts,+ FlexibleInstances,+ UndecidableInstances+ other-extensions:+ TemplateHaskell+ ghc-options: -Wall+ build-depends:+ aern2-fun >=0.2.9+ , aern2-mp >=0.2.9+ , base ==4.*+ , collect-errors >=0.1.5+ , mixed-types-num >=0.5.9+ , vector+ default-language: Haskell2010++executable aern2-mfun-benchmark-optimisation+ main-is: bench/OptimisationBenchmark.hs+ other-modules:+ Paths_aern2_mfun+ default-extensions:+ RebindableSyntax,+ ScopedTypeVariables,+ TypeFamilies,+ TypeOperators,+ ConstraintKinds,+ DefaultSignatures,+ MultiParamTypeClasses,+ FlexibleContexts,+ FlexibleInstances,+ UndecidableInstances+ other-extensions:+ TemplateHaskell+ ghc-options: -threaded -rtsopts -with-rtsopts=-N -Wall -O2+ build-depends:+ aern2-fun >=0.2.9+ , aern2-mfun+ , aern2-mp >=0.2.9+ , base ==4.*+ , collect-errors >=0.1.5+ , mixed-types-num >=0.5.9+ , vector+ default-language: Haskell2010
+ bench/OptimisationBenchmark.hs view
@@ -0,0 +1,30 @@+module Main where++import MixedTypesNumPrelude+import AERN2.BoxFun.TestFunctions+import AERN2.BoxFun.Optimisation+import AERN2.MP.Ball+-- import AERN2.Real+-- import AERN2.MP.Dyadic+-- import AERN2.MP.Accuracy+-- import AERN2.AD.Differential++-- import qualified Data.List as List++-- import AERN2.BoxFun.Type+-- import AERN2.Linear.Vector.Type as V++main :: IO ()+main = + let+ ac = bits 100+ in+ do+ putStrLn $ "ratz: " ++ (show $ minFun ratz4 ac (prec 53))+ putStrLn $ "griewank 2: " ++ (show $ minFun (griewank 2) ac (prec 53))+ putStrLn $ "shekel 2: " ++ (show $ minFun shekel ac (prec 53))+ putStrLn $ "himmelblau: " ++ (show $ minFun himmelblau ac (prec 53))+ putStrLn $ "rosenbrock: " ++ (show $ minFun rosenbrock ac (prec 53))+ putStrLn $ "trefethen: " ++ (show $ minFun siam4 ac (prec 53))+ putStrLn $ "griewank 5: " ++ (show $ minFun (griewank 5) ac (prec 53))+ putStrLn $ "griewank 7: " ++ (show $ minFun (griewank 7) ac (prec 53))
+ src/AERN2/AD/Differential.hs view
@@ -0,0 +1,11 @@+module AERN2.AD.Differential +(+ module AERN2.AD.Type+ , module AERN2.AD.GenericOperations+ , module AERN2.AD.MPBallOperations+)+where++import AERN2.AD.Type+import AERN2.AD.GenericOperations+import AERN2.AD.MPBallOperations
+ src/AERN2/AD/GenericOperations.hs view
@@ -0,0 +1,220 @@+{-# OPTIONS_GHC -Wno-orphans #-}+module AERN2.AD.GenericOperations where++import MixedTypesNumPrelude+import AERN2.AD.Type++instance + (CanDiv a a, CanSubSameType a, CanMulSameType a, CanMulBy a Integer, CanAddSameType a,+ CanSubSameType (DivType a a)) + =>+ CanDiv (Differential a) (Differential a)+ where+ type DivType (Differential a) (Differential a) = Differential (DivType a a)+ divide a b =+ case min (order a) (order b) of+ 2 -> + let+ dtDiff = (a_dxt * b_x - a_x * b_dxt)+ ySqrd = b_x * b_x+ exp0 = a_d2x * b_x + a_dxt * b_dx - b_dxt * a_dx - a_x * b_d2x+ in+ OrderTwo (a_x / b_x) ((a_dx * b_x - a_x * b_dx)/ySqrd) (dtDiff/ySqrd) + (exp0/ySqrd - (2*b_dx * dtDiff) / (ySqrd * b_x))+ 1 -> OrderOne (a_x / b_x) ((a_dx * b_x - a_x * b_dx)/(b_x * b_x)) + 0 -> OrderZero (a_x / b_x)+ _ -> error "illegal Differential order"+ where+ a_x = diff_x a+ b_x = diff_x b+ a_dx = diff_dx a+ b_dx = diff_dx b+ a_dxt = diff_dxt a+ b_dxt = diff_dxt b+ a_d2x = diff_d2x a+ b_d2x = diff_d2x b++instance + (CanExpSameType a, CanMulSameType a, CanAddSameType a) =>+ CanExp (Differential a)+ where+ type ExpType (Differential a) = Differential a+ exp (OrderZero x) = OrderZero (exp x)+ exp (OrderOne x dx) = OrderOne (exp x) (dx * exp x)+ exp (OrderTwo x dx dxt d2x) = OrderTwo (exp x) (dx * exp x) (dxt * exp x) ((dxt * dx + d2x) * exp x)+++{- TODO: fix this in AERN2.MP -}+clampedCos :: (CanSinCosSameType a, CanMinMaxSameType a, HasIntegers a) => a -> a+clampedCos (x :: a) = max (convertExactly $ -1 :: a) $ min ((convertExactly 1) :: a) (cos x)++clampedSin :: (CanSinCosSameType a, CanMinMaxSameType a, HasIntegers a) => a -> a+clampedSin (x :: a) = max (convertExactly $ -1 :: a) $ min ((convertExactly 1) :: a) (sin x)++instance + (CanSinCosSameType a, CanMulSameType a, CanNegSameType a, CanSubSameType a, CanAddSameType a+ , HasIntegers a, CanMinMaxSameType a) + =>+ CanSinCos (Differential a)+ where+ type SinCosType (Differential a) = Differential a+ cos (OrderZero x) = OrderZero (clampedCos x)+ cos (OrderOne x dx) = OrderOne (clampedCos x) (-dx * clampedSin x)+ cos (OrderTwo x dx dxt d2x) = OrderTwo (clampedCos x) (-dx * clampedSin x) (-dxt * clampedSin x) + (-dxt * dx * clampedCos x - d2x * clampedSin x)+ + sin (OrderZero x) = OrderZero (clampedSin x)+ sin (OrderOne x dx) = OrderOne (clampedSin x) (dx * clampedCos x)+ sin (OrderTwo x dx dxt d2x) = OrderTwo (clampedSin x) (dx * clampedCos x) (dxt * clampedCos x)+ (d2x * clampedCos x - dxt * dx * clampedSin x)++instance + (CanAddSameType a) => + CanAddAsymmetric (Differential a) (Differential a)+ where+ type AddType (Differential a) (Differential a) = Differential a+ add a b =+ case min (order a) (order b) of+ 2 -> OrderTwo (a_x + b_x) (a_dx + b_dx) (a_dxt + b_dxt) (a_d2x + b_d2x)+ 1 -> OrderOne (a_x + b_x) (a_dx + b_dx)+ 0 -> OrderZero (a_x + b_x)+ _ -> error "illegal Differential order"+ where+ a_x = diff_x a+ b_x = diff_x b+ a_dx = diff_dx a+ b_dx = diff_dx b+ a_dxt = diff_dxt a+ b_dxt = diff_dxt b+ a_d2x = diff_d2x a+ b_d2x = diff_d2x b+ +instance + (CanMulSameType a, CanAddSameType a) =>+ CanMulAsymmetric (Differential a) (Differential a)+ where+ type MulType (Differential a) (Differential a) = Differential a+ mul a b =+ case min (order a) (order b) of+ 2 -> OrderTwo (a_x * b_x) (a_dx * b_x + a_x * b_dx) (a_x * b_dxt + b_x * a_dxt) + (a_x * b_d2x + a_dx * b_dxt + a_dxt * b_dx + a_d2x * b_x)+ 1 -> OrderOne (a_x * b_x) (a_dx * b_x + a_x * b_dx)+ 0 -> OrderZero (a_x * b_x)+ _ -> error "illegal Differential order"+ where+ a_x = diff_x a+ b_x = diff_x b+ a_dx = diff_dx a+ b_dx = diff_dx b+ a_dxt = diff_dxt a+ b_dxt = diff_dxt b+ a_d2x = diff_d2x a+ b_d2x = diff_d2x b++instance + (CanMulSameType a, CanAddSameType a, CanPowBy a a, CanSubThis a Integer, CanLogSameType a, CanDivSameType a) =>+ CanPow (Differential a) (Differential a)+ where+ type PowType (Differential a) (Differential a) = (Differential a)+ pow a b =+ case min (order a) (order b) of+ 2 -> OrderOne (a_x ^ b_x) ((a_x ^ (b_x - 1)) * (b_x * a_dx + a_x * (log a_x) * b_dx)) --FIXME: Add real OrderTwo definition here+ -- 2 -> OrderTwo (x a ^ x b) + -- ((x a ^ (x b - (convertExactly 1 :: a))) * (x b * dx a + x a * log (x a) * dx b)) + -- ((x a ^ (x b - (convertExactly 1 :: a))) * (x b * dxt a + x a * log (x a) * dxt b))+ -- (x a)+ + -- where+ -- ta = convertExactly 2 :: a+ 1 -> OrderOne (a_x ^ b_x) ((a_x ^ (b_x - 1)) * (b_x * a_dx + a_x * (log a_x) * b_dx))+ 0 -> OrderZero (a_x ^ b_x)+ _ -> undefined+ where+ a_x = diff_x a+ b_x = diff_x b+ a_dx = diff_dx a+ b_dx = diff_dx b+ -- a_dxt = diff_dxt a+ -- b_dxt = diff_dxt b+ -- a_d2x = diff_d2x a+ -- b_d2x = diff_d2x b++instance + (CanSubSameType a) =>+ CanSub (Differential a) (Differential a)+ where+ type SubType (Differential a) (Differential a) = Differential a+ sub a b =+ case min (order a) (order b) of+ 2 -> OrderTwo (a_x - b_x) (a_dx - b_dx) (a_dxt - b_dxt) (a_d2x - b_d2x)+ 1 -> OrderOne (a_x - b_x) (a_dx - b_dx)+ 0 -> OrderZero (a_x - b_x)+ _ -> error "illegal Differential order"+ where+ a_x = diff_x a+ b_x = diff_x b+ a_dx = diff_dx a+ b_dx = diff_dx b+ a_dxt = diff_dxt a+ b_dxt = diff_dxt b+ a_d2x = diff_d2x a+ b_d2x = diff_d2x b++instance + (CanMulBy (Differential a) Integer) =>+ CanNeg (Differential a)+ where+ type NegType (Differential a) = Differential a+ negate x = (-1) * x++instance + (CanSqrtSameType a, CanMulSameType a, CanNegSameType a, CanAddSameType a, CanMulBy a Integer, CanRecipSameType a) + =>+ CanSqrt (Differential a)+ where+ type SqrtType (Differential a) = Differential a+ sqrt (OrderZero x) = OrderZero (sqrt x)+ sqrt (OrderOne x dx) = OrderOne (sqrt x) (dx * sqrtx')+ where sqrtx' = recip (2 * sqrt x)+ sqrt (OrderTwo x dx dxt d2x) = OrderTwo (sqrt x) (dx * sqrtx') (dxt * sqrtx') + ((d2x * sqrtx') + (dx * dxt * sqrtx'')) + where + sqrtx' = recip (2 * sqrt x)+ sqrtx'' = negate $ recip (4 * x * sqrt x)+ -- sqrtx' == 1 / (2 * sqrt(x))+ -- sqrtx'' == -1 / (4 * x * sqrt(x)) == -1 / (4 * x^(3/2))++instance+ CanMinMaxSameType a =>+ CanMinMaxAsymmetric (Differential a) (Differential a)+ where+ type MinMaxType (Differential a) (Differential a) = Differential a+ min a b = OrderZero (min (diff_x a) (diff_x b))+ max a b = OrderZero (max (diff_x a) (diff_x b))+++-- instance+-- (CanMinMaxSameType a, HasIntegers a) =>+-- CanMinMaxAsymmetric (Differential a) (Differential a)+-- where+-- type MinMaxType (Differential a) (Differential a) = Differential a+-- min a b = +-- case min (order a) (order b) of+-- 2 -> OrderTwo (min (a_x) (b_x)) (min (a_dx) (b_dx)) (min (a_dxt) (b_dxt)) +-- (min (a_d2x) (b_d2x))+-- 1 -> OrderOne (min (a_x) (b_x)) (min (a_dx) (b_dx))+-- 0 -> OrderZero (min (a_x) (b_x))+-- max a b =+-- case min (order a) (order b) of+-- 2 -> OrderTwo (max (a_x) (b_x)) (max (a_dx) (b_dx)) (max (a_dxt) (b_dxt)) +-- (max (a_d2x) (b_d2x))+-- 1 -> OrderOne (max (a_x) (b_x)) (max (a_dx) (b_dx))+-- 0 -> OrderZero (max (a_x) (b_x))++-- instance+-- (CanAbsSameType a) =>+-- CanAbs (Differential a)+-- where+-- type AbsType (Differential a) = Differential a+-- abs = fmap abs+
+ src/AERN2/AD/MPBallOperations.hs view
@@ -0,0 +1,223 @@+{-# OPTIONS_GHC -Wno-orphans #-}+module AERN2.AD.MPBallOperations where++import MixedTypesNumPrelude++import AERN2.AD.Type+import AERN2.AD.GenericOperations ()+import AERN2.MP.Ball+import Numeric.CollectErrors.Type (noValueNumErrorCertain, NumError (NumError))++{-- addition --}++instance+ CanAddAsymmetric MPBall (Differential (CN MPBall))+ where+ type AddType MPBall (Differential (CN MPBall)) = Differential (CN MPBall)+ add a b = add (differential (order b) (cn a)) b++instance + CanAddAsymmetric (Differential (CN MPBall)) MPBall+ where+ type AddType (Differential (CN MPBall)) MPBall = Differential (CN MPBall)+ add b a = add b (differential (order b) (cn a))++instance + CanAddAsymmetric Integer (Differential (CN MPBall))+ where+ type AddType Integer (Differential (CN MPBall)) = Differential (CN MPBall)+ add a b = add (differential (order b) (cn $ mpBall a)) b++instance + CanAddAsymmetric (Differential (CN MPBall)) Integer+ where+ type AddType (Differential (CN MPBall)) Integer = Differential (CN MPBall)+ add b a = add b (differential (order b) (cn $ mpBall a))++instance+ (CanBeMPBall a) =>+ CanAddAsymmetric (CN a) (Differential (CN MPBall))+ where+ type AddType (CN a) (Differential (CN MPBall)) = Differential (CN MPBall)+ add a b = add (differential (order b) (fmap mpBall a)) b++instance + (CanBeMPBall a) =>+ CanAddAsymmetric (Differential (CN MPBall)) (CN a)+ where+ type AddType (Differential (CN MPBall)) (CN a) = Differential (CN MPBall)+ add b a = add b (differential (order b) (fmap mpBall a))++{-- subtraction --}++instance + CanSub Integer (Differential (CN MPBall))+ where+ type SubType Integer (Differential (CN MPBall)) = Differential (CN MPBall)+ sub a b = sub (differential (order b) (cn $ mpBall a)) b++instance + CanSub (Differential (CN MPBall)) Integer+ where+ type SubType (Differential (CN MPBall)) Integer = Differential (CN MPBall)+ sub b a = sub b (differential (order b) (cn $ mpBall a))++instance + CanSub MPBall (Differential (CN MPBall))+ where+ type SubType MPBall (Differential (CN MPBall)) = Differential (CN MPBall)+ sub a b = sub (differential (order b) (cn a)) b++instance + CanSub (Differential (CN MPBall)) MPBall+ where+ type SubType (Differential (CN MPBall)) MPBall = Differential (CN MPBall)+ sub b a = sub b (differential (order b) (cn a))++instance + (CanBeMPBall a) =>+ CanSub (CN a) (Differential (CN MPBall))+ where+ type SubType (CN a) (Differential (CN MPBall)) = Differential (CN MPBall)+ sub a b = sub (differential (order b) (fmap mpBall a)) b++instance + (CanBeMPBall a) =>+ CanSub (Differential (CN MPBall)) (CN a) + where+ type SubType (Differential (CN MPBall)) (CN a) = Differential (CN MPBall)+ sub b a = sub b (differential (order b) (fmap mpBall a))++{-- multiplication --}++instance + CanMulAsymmetric Integer (Differential (CN MPBall))+ where+ type MulType Integer (Differential (CN MPBall)) = Differential (CN MPBall)+ mul a b = mul (differential (order b) (cn $ mpBall a)) b++instance + CanMulAsymmetric (Differential (CN MPBall)) Integer+ where+ type MulType (Differential (CN MPBall)) Integer = Differential (CN MPBall)+ mul b a = mul b (differential (order b) (cn $ mpBall a))++instance + CanMulAsymmetric MPBall (Differential (CN MPBall))+ where+ type MulType MPBall (Differential (CN MPBall)) = Differential (CN MPBall)+ mul a b = mul (differential (order b) (cn a)) b++instance + CanMulAsymmetric (Differential (CN MPBall)) MPBall+ where+ type MulType (Differential (CN MPBall)) MPBall = Differential (CN MPBall)+ mul b a = mul b (differential (order b) (cn a))++{-- division --}++instance + CanDiv Integer (Differential (CN MPBall))+ where+ type DivType Integer (Differential (CN MPBall)) = Differential (CN MPBall)+ divide a b = divide (differential (order b) (cn $ mpBall a)) b+ +instance + CanDiv (Differential (CN MPBall)) Integer+ where+ type DivType (Differential (CN MPBall)) Integer = Differential (CN MPBall)+ divide b a = divide b (differential (order b) (cn $ mpBall a)) ++instance + CanDiv MPBall (Differential (CN MPBall))+ where+ type DivType MPBall (Differential (CN MPBall)) = Differential (CN MPBall)+ divide a b = divide (differential (order b) (cn a)) b+ +instance + CanDiv (Differential (CN MPBall)) MPBall+ where+ type DivType (Differential (CN MPBall)) MPBall = Differential (CN MPBall)+ divide b a = divide b (differential (order b) (cn a)) ++instance + (CanBeMPBall a) =>+ CanMulAsymmetric (CN a) (Differential (CN MPBall))+ where+ type MulType (CN a) (Differential (CN MPBall)) = Differential (CN MPBall)+ mul a b = mul (differential (order b) (fmap mpBall a)) b++instance + (CanBeMPBall a) =>+ CanMulAsymmetric (Differential (CN MPBall)) (CN a)+ where+ type MulType (Differential (CN MPBall)) (CN a) = Differential (CN MPBall)+ mul b a = mul b (differential (order b) (fmap mpBall a))++instance + (CanBeMPBall a) =>+ CanDiv (CN a) (Differential (CN MPBall))+ where+ type DivType (CN a) (Differential (CN MPBall)) = Differential (CN MPBall)+ divide a b = divide (differential (order b) (fmap mpBall a)) b++instance + (CanBeMPBall a) =>+ CanDiv (Differential (CN MPBall)) (CN a) + where+ type DivType (Differential (CN MPBall)) (CN a) = Differential (CN MPBall)+ divide b a = divide b (differential (order b) (fmap mpBall a))+++instance + CanPow (Differential (CN MPBall)) Integer+ where+ type PowType (Differential (CN MPBall)) Integer = (Differential (CN MPBall))+ pow (OrderZero x) n =+ OrderZero p+ where+ rawP = pow x n+ p = if even n then abs rawP else rawP -- TODO: this should go in MPBall+ pow (OrderOne x dx) n = + OrderOne p dp+ where+ nEven = even n+ rawP' = pow x (n - 1)+ rawP = x * rawP'+ p = if nEven then abs rawP else rawP -- TODO: this should go in MPBall+ dp = n * dx * if nEven then rawP' else abs rawP'+ pow (OrderTwo x dx dxt d2x) n =+ OrderTwo p dp dpt d2p+ where+ nEven = even n+ rawP'' = pow x (n - 2)+ rawP' = x * rawP''+ rawP = x * rawP'+ p = if nEven then abs rawP else rawP -- TODO: this should go in MPBall+ xnm1 = if nEven then rawP' else abs rawP'+ dp = n * dx * xnm1+ dpt = n * dxt * xnm1+ xnm2 = if nEven then abs rawP'' else rawP''+ d2p = n*((n - 1)*dx*dxt*xnm2 + d2x*xnm1)++instance+ CanAbs (Differential (CN MPBall))+ where+ type AbsType (Differential (CN MPBall)) = Differential (CN MPBall)+ abs (OrderZero x) = OrderZero $ abs x + abs (OrderOne x dx) = OrderOne (abs x) newDx+ where+ newDx = do + dx_ <- dx+ pure $ (hullMPBall dx_ (-dx_))+ abs (OrderTwo _ _ _ _) = error "Abs for differential order two undefined"++instance + CanDivIMod (Differential (CN MPBall)) (Differential (CN MPBall))+ where+ type DivIType (Differential (CN MPBall)) (Differential (CN MPBall)) = (Differential (CN MPBall))+ divIMod a b = (error "Integer division for Differential (CN MPBall) undefined", OrderTwo (mod ax bx) err err err)+ where + ax = diff_x a+ bx = diff_x b+ err = noValueNumErrorCertain $ NumError "No derivatives after modulus"
+ src/AERN2/AD/Type.hs view
@@ -0,0 +1,46 @@+module AERN2.AD.Type where++import MixedTypesNumPrelude+import AERN2.MP.Precision++data Differential a =+ OrderZero {diff_x :: a}+ | OrderOne {diff_x :: a, diff_dx :: a}+ | OrderTwo {diff_x :: a, diff_dx :: a, diff_dxt :: a, diff_d2x :: a}+ deriving (Show)++order :: Differential a -> Integer+order (OrderZero _) = 0+order (OrderOne _ _) = 1+order (OrderTwo _ _ _ _) = 2++class CanBeDifferential a where+ differential :: Integer -> a -> Differential a++instance + (HasIntegers a) =>+ CanBeDifferential a+ where + differential 0 a = OrderZero a+ differential 1 a = OrderOne a (convertExactly 0)+ differential _ a = OrderTwo a (convertExactly 0) (convertExactly 0) (convertExactly 0)++instance Functor Differential where+ fmap f (OrderZero x) = OrderZero (f x)+ fmap f (OrderOne x dx) = OrderOne (f x) (f dx)+ fmap f (OrderTwo x dx dxt d2x) = OrderTwo (f x) (f dx) (f dxt) (f d2x)++instance + (HasPrecision a) => (HasPrecision (Differential a))+ where+ getPrecision a = getPrecision (diff_x a) -- TODO: safe?++instance + (CanSetPrecision a) => (CanSetPrecision (Differential a))+ where+ setPrecision p = fmap (setPrecision p)++setValue :: Differential a -> a -> Differential a+setValue (OrderZero _x) v = OrderZero v+setValue (OrderOne _x dx) v = OrderOne v dx+setValue (OrderTwo _x dx dxt d2x) v = OrderTwo v dx dxt d2x
+ src/AERN2/BoxFun/Box.hs view
@@ -0,0 +1,150 @@+{-# OPTIONS_GHC -Wno-orphans #-}+module AERN2.BoxFun.Box where++import qualified Prelude as Prelude+import MixedTypesNumPrelude+import AERN2.MP.Dyadic+import AERN2.MP.Ball+import AERN2.Linear.Vector.Type (Vector, (!))+import qualified AERN2.Linear.Vector.Type as V+import AERN2.Util.Util++type Box = Vector (CN MPBall)++-- instance HasEqAsymmetric Box Box where+-- type EqCompareType Box Box = (CN Kleenean)+-- equalTo box0 box1 = +-- V.foldl' (&&) (cn CertainTrue) $ V.zipWith equalTo box0 box1++-- intersection :: Box -> Box -> Maybe Box+-- intersection = undefined++intersectionCertainlyEmpty :: Box -> Box -> Bool+intersectionCertainlyEmpty vx vy =+ V.foldl' (||) False $ V.zipWith ballIntersectionCertainlyEmpty vx vy+ where + ballIntersectionCertainlyEmpty x y = + (lx !>! ry || rx !<! ly)+ where+ (lx, rx) = endpointsAsIntervals x+ (ly, ry) = endpointsAsIntervals y++nonEmptyIntersection :: Box -> Box -> Box+nonEmptyIntersection vx vy = + V.zipWith intersectCN vx vy++instance IsBall Box where+ type CentreType (Box) = Box+ centre v = V.map (\x -> (cnMPBallP (getPrecision x) . centre) x) v++inftyNorm :: Box -> CN MPBall+inftyNorm box = + V.foldl' (\n x -> max n (abs x)) (cn $ mpBall 0) box++ellOneNorm :: Box -> CN MPBall+ellOneNorm box = + V.foldl' (\n x -> n + (abs x)) (cn $ mpBall 0) box++width :: Box -> CN MPBall+width box = + V.foldl' (\n x -> max n (2 * ((fmap mpBall) $ (fmap radius) x) ) ) (cn $ mpBall 0) box++widestDirection :: Box -> Integer+widestDirection box =+ aux (V.length box - 1) (V.length box - 1) (dyadic 0)+ where+ aux k i lth =+ if k < 0 then + i+ else+ let + x = box ! k+ lth' = + 2 * (dyadic $ radius $ unCN x) -- TODO: unsafe+ in+ if (lth' > lth) then + aux (k - 1) k lth'+ else + aux (k - 1) i lth+++bisect :: Integer -> Box -> (Box, Box)+bisect k box =+ if exponent_ > 0+ && (abs $ mc - lbc) < (dyadic 0.5)^exponent_ then+ (setPrecision (increasePrecision p) lb, setPrecision (increasePrecision p) rb)+ else + (lb, rb)+ where+ exponent_ = integer p - ilog+ NormBits ilog = 2 + (getNormLog $ mpBall $ 1 + abs lbc )+ lbc = unCN (centre $ lb ! k)+ -- rbc = unCN (centre $ rb ! k)+ mc = unCN (centre $ m)+ increasePrecision p1 =+ p1 + (prec $ (integer p1) `Prelude.div` 2)+ lb = setPrecision p leftBox+ rb = setPrecision p rightBox+ p = getPrecision box+ interval = box ! k+ m = (fmap $ mpBallP p) $ (fmap centre) interval+ (l, r) = endpointsAsIntervals interval+ leftBox = V.imap (\i x -> if i == k then fromEndpointsAsIntervals l m else x) box+ rightBox = V.imap (\i x -> if i == k then fromEndpointsAsIntervals m r else x) box++-- Bisects a into 2^d boxes of the same size, where d is the dimension of the given box+fullBisect :: Box -> [Box]+fullBisect b =+ case V.length b of+ 0 -> [b]+ l ->+ -- y is the dimension bisected in the current iteration+ -- x is a bisection of the previous dimension (tail recursion)+ concatMap (\x -> map (\y -> x V.+++ V.singleton y) (bisectDimension (l-1))) (fullBisect (V.take (fromIntegral (l-1)) b))++ where+ bisectDimension n = [fst bn ! n, snd bn ! n]+ where bn = bisect n b++-- Get the endpoints of a box as a list containing a pair of MPBalls for each dimension+getEndpoints :: Box -> [(MPBall, MPBall)]+getEndpoints b = + case V.length b of+ 0 -> []+ _ -> endpointsAsIntervals (unCN (b ! 0)) : getEndpoints (V.drop (int 1) b)++lowerBounds :: Box -> Box+lowerBounds = V.map lowerBound++upperBounds :: Box -> Box+upperBounds = V.map upperBound++instance + CanNeg Box+ where+ type NegType Box = Box+ negate box = V.map (\x -> -x) box+ +createEnclosingBox :: Box -> Box -> Box+createEnclosingBox box1 box2 =+ enclosingBox+ where+ indexedBox1 = V.zip (V.fromList [0 .. integer (V.length box1) - 1]) box1+ enclosingBox = + V.map + (\(i, x) -> + let+ y = box2 V.! i+ (l, r) = endpointsAsIntervals x+ (l', r') = endpointsAsIntervals y+ newL = min l l'+ newR = max r r'+ in+ fromEndpointsAsIntervals newL newR)+ indexedBox1++intersectList :: [Box] -> Box+intersectList [] = V.empty+intersectList [b] = b+intersectList [b1,b2] = if intersectionCertainlyEmpty b1 b2 then V.empty else nonEmptyIntersection b1 b2+intersectList (b1:b2:bs) = if intersectionCertainlyEmpty b1 b2 then V.empty else intersectList $ nonEmptyIntersection b1 b2 : bs
+ src/AERN2/BoxFun/Optimisation.hs view
@@ -0,0 +1,319 @@+module AERN2.BoxFun.Optimisation where++import qualified Prelude+import MixedTypesNumPrelude+import qualified Numeric.CollectErrors as CN+import AERN2.MP.Dyadic+import AERN2.MP.Ball+import AERN2.BoxFun.Box (Box)+import qualified AERN2.BoxFun.Box as Box+import AERN2.BoxFun.Type+import AERN2.Kleenean+import AERN2.Linear.Vector.Type as V+import AERN2.Linear.Matrix.Type+import AERN2.Linear.Matrix.Inverse++import qualified Data.List as List++import qualified AERN2.PQueue as Q++import AERN2.Util.Util++import Debug.Trace (trace)++globalMinimumGreaterThanN :: BoxFun -> Accuracy -> CN Rational -> Precision -> Bool+globalMinimumGreaterThanN f ac n initialPrecision =+ trace (show x)+ x !>! n+ where x = globalMinimum f ac initialPrecision ++minFun :: BoxFun -> Accuracy -> Precision -> (Integer, CN MPBall)+minFun f ac initialPrecision = + bestLocalMinimum f (domain f) ac initialPrecision++data SearchBox =+ SearchBox+ {+ extents :: Box+ , minimum :: CN MPBall + } deriving (Show)++instance + HasPrecision SearchBox+ where+ getPrecision (SearchBox b _) = getPrecision b++instance + CanSetPrecision SearchBox+ where+ setPrecision p (SearchBox b m) = SearchBox (setPrecision p b) m++instance Prelude.Eq SearchBox where+ (==) (SearchBox _ _) (SearchBox _ _) =+ False -- TODO: safe?++instance Prelude.Ord SearchBox where+ (<=) (SearchBox _ min0) (SearchBox _ min1) = + case (CN.toEither $ (lowerBound min0 :: CN MPBall), CN.toEither $ (lowerBound min1 :: CN MPBall)) of+ (Left _, Left _) -> True+ (Left _, Right _ ) -> True+ (Right _ , Left _) -> False+ (Right m0, Right m1) -> + centre m0 - (dyadic $ radius m0) <= centre m1 - (dyadic $ radius m1) -- TODO: radius should be 0++---++globalMinimumWithCutoff :: BoxFun -> Accuracy -> CN MPBall -> Precision -> CN MPBall+globalMinimumWithCutoff f ac cutoff initialPrecision =+ if dimension f == 1 then+ let+ fl = apply f (V.map lowerBound $ domain f)+ fr = apply f (V.map upperBound $ domain f)+ localMin = snd $ bestLocalMinimumWithCutoff f (domain f) ac cutoff initialPrecision+ in+ min fl $ min localMin fr+ else + let+ localMin = snd $ bestLocalMinimumWithCutoff f (domain f) ac cutoff initialPrecision+ boundaryFuns = boundaryRestrictions f+ boundaryMinima = List.map (\g -> globalMinimumWithCutoff g ac (min cutoff ((upperBound localMin :: CN MPBall))) initialPrecision) boundaryFuns+ in+ List.foldl' min localMin boundaryMinima+++globalMinimum :: BoxFun -> Accuracy -> Precision -> CN MPBall+globalMinimum f ac initialPrecision =+ globalMinimumWithCutoff f ac (apply f (centre boxp)) initialPrecision+ where+ boxp = setPrecision initialPrecision (domain f)++bestLocalMinimum :: BoxFun -> Box -> Accuracy -> Precision -> (Integer, CN MPBall)+bestLocalMinimum f box ac initialPrecision =+ bestLocalMinimumWithCutoff f box ac (apply f (centre boxp)) initialPrecision+ where+ boxp = setPrecision initialPrecision box++bestLocalMinimumWithCutoff :: BoxFun -> Box -> Accuracy -> CN MPBall -> Precision -> (Integer, CN MPBall)+bestLocalMinimumWithCutoff f box ac initialCutoff initialPrecision =+ aux initialQueue initialCutoff 0 dummyBox+ where+ boxp = setPrecision initialPrecision box+ initialRange = apply f boxp+ initialSearchBox = SearchBox boxp initialRange+ initialQueue = Q.singleton initialSearchBox+ dummyBox = SearchBox (V.fromList [cn $ mpBall $ 10^6]) initialRange -- TODO: hack...++ aux q cutoff steps (SearchBox _lastBox rng) = + case Q.minView q of+ Nothing -> trace ("no local minimum.") $ (steps, rng)+ Just (minBox, q') ->+ --trace ("value: "++ (show $ val)) $+ trace ("min box: "++ (show $ minBox)) $+ --trace ("box acc: "++ (show $ getAccuracy $ ext)) $+ --trace (show $ Box.width (extents minBox)) $+ --trace ("lower bound "++ (show $ Box.lowerBound $ val)) $+ --trace ("val' "++ (show $ val')) $+ trace ("cutoff: "++ (show $ cutoff)) $+ trace ("queue size: "++ (show $ Q.size q)) $+ --trace ("cutoff == 0? "++(show $ cutoff == (mpBall 0))) $+ --trace ("precision: "++ (show $ precision)) $+ --trace ("dist to last "++ (show $ distToLast)) $+ --trace ("accuracy: "++ (show $ getAccuracy val')) $+ --trace ("precision centre: "++ (show $ fmap (getPrecision . centre) val)) $+ if getAccuracy val' >= ac then+ (steps, val')+ else + aux q'' newCutoff (steps + 1) (SearchBox ext rng)+ where+ val' = fromEndpointsAsIntervals (lowerBound val) (cutoff)+ SearchBox ext val = minBox++ (newCutoff, newBoxes) = + processBox f ac cutoff minBox++ q'' = foldr (Q.insert) q' newBoxes++lipschitzContraction :: BoxFun -> Box -> SearchBox -> SearchBox+lipschitzContraction f g (SearchBox box m) =+ {-trace("fa: "++(show $ getAccuracy (apply f box))) $+ trace("la: "++(show $ getAccuracy $ dotProduct)) $+ trace("ba: "++(show $ getAccuracy $ box ! int 0)) $-}+ {-if (radius $ (~!) $ newRange) < (radius $ (~!) $ m) then+ trace ("Lipschitz better.")+ SearchBox box m'+ else -}+ SearchBox box m'+ where+ boxCentre = centre box+ centreValue = apply f boxCentre+ difference = box - boxCentre+ dotProduct = g * difference + newRange = centreValue + dotProduct+ m' = intersectCN m newRange++lipschitzRange :: BoxFun -> CN MPBall -> Box -> Box -> Box -> CN MPBall -> CN MPBall+lipschitzRange _f fc c g box m =+ m'+ where+ difference = box - c+ normG = Box.ellOneNorm g+ normDiff = Box.inftyNorm difference+ dotProduct = normG * normDiff+ newRange = fc + (fromEndpointsAsIntervals (-dotProduct) dotProduct :: CN MPBall)+ m' = intersectCN m newRange++applyLipschitz :: BoxFun -> Box -> CN MPBall+applyLipschitz f box = + lipschitzRange f fbc bc dfb' box fb+ where+ (fb, dfb') = valueGradient f box+ bc = centre box+ fbc = apply f bc++increasePrecision :: Precision -> Precision+increasePrecision p =+ p + (prec $ (integer p) `Prelude.div` 2)++newtonStep :: BoxFun -> Accuracy -> Vector (CN MPBall) -> Vector (CN MPBall) -> Matrix (CN MPBall) -> SearchBox -> Bool -> Maybe (Bool, SearchBox)+newtonStep f ac c dfc hInv b@(SearchBox box m) newtonSuccesful = + --Just $ SearchBox box' m'+ {-trace ("precision m "++(show $ (fmap getPrecision) m)) $+ trace ("precision m' "++(show $ (fmap getPrecision) m')) $+ trace ("precision box centre "++(show $ getPrecision c)) $+ trace ("precision box "++(show $ getPrecision box)) $+ trace ("precision newton box "++(show $ getPrecision newtonBox)) $+ trace ("precision box' "++(show $ getPrecision box')) $+ trace ("precision hInv "++(show $ getPrecision (entries hInv ! int 0))) $-}+ if getAccuracy m >= ac then+ Just (newtonSuccesful, b)+ --else if not hInvDefined then+ -- Just (newtonSuccesful, b)+ else if Box.intersectionCertainlyEmpty box newtonBox then+ Nothing+ else if Box.width box' !<=! (dyadic $ 0.75) * Box.width box then+ if getAccuracy m' > getAccuracy m then+ newtonStep f ac c dfc hInv (SearchBox box' m') True+ else + Just (True, SearchBox (setPrecision (increasePrecision $ getPrecision box') box') m')+ else + Just (newtonSuccesful, SearchBox box' m')+ where+ {-c = centre box+ dfc = gradient f c-}+ -- hInvDefined = V.foldl' (&&) (True) $ V.map (isJust . fst . ensureNoCN) (entries hInv)+ newtonBox = c - hInv * (dfc)+ box' = Box.nonEmptyIntersection box newtonBox+ m' = apply f box'++processBox :: BoxFun -> Accuracy -> CN MPBall -> SearchBox -> (CN MPBall, [SearchBox])+processBox f ac cutoff box =+ if getAccuracy ext < bits 10 then + split f (gradient f ext) cutoff ext + else + result+ where+ ext = extents box+ (_fb, dfb, hfb) = valueGradientHessian f ext+ c = centre ext+ dfc = gradient f c+ maybeHinv = inverse hfb+ -- p = getPrecision box+ box' = --Just (False, box)+ case maybeHinv of + Nothing -> Just (False, box)+ Just hInv -> newtonStep f ac c dfc hInv box False+ result =+ case box' of + Nothing -> (cutoff, [])+ Just (newtonSuccesful, bx@(SearchBox bxe m)) ->+ let+ c' = min (upperBound $ apply f $ centre bxe :: CN MPBall) cutoff + in+ if newtonSuccesful then+ if getAccuracy m >= ac then+ (c', [bx])+ else+ processBox f ac c' bx+ else+ split f dfb c' bxe++split :: BoxFun -> Vector (CN MPBall) -> CN MPBall -> Box -> (CN MPBall, [SearchBox])+split f dfb cutoff bxe = + let+ diff = bxe - centre bxe+ dir i = (fmap dyadic) $ (fmap radius) $ (dfb ! i) * (diff ! i) :: CN Dyadic+ dirs = V.map dir $ V.enumFromTo 0 (V.length bxe - 1)+ dirsDefined = V.foldl' (&&) True $ V.map (not . CN.hasError) dirs+ aux k j d = + if k == V.length bxe then + j + else + let+ d' = unCN $ dirs ! k+ in+ if d' > d then + aux (k + 1) k d'+ else+ aux (k + 1) j d+ splittingIndex = + if dirsDefined then (aux 1 0 (unCN $ dirs ! 0)) else Box.widestDirection bxe+ (a , b) = Box.bisect splittingIndex bxe+ (fa, dfa') = valueGradient f a+ (fb, dfb') = valueGradient f b+ ac = centre a+ bc = centre b+ fac = apply f ac+ fbc = apply f bc+ fa' = lipschitzRange f fac ac dfa' a fa+ fb' = lipschitzRange f fbc bc dfb' b fb+ cutoff' = min (upperBound fac :: CN MPBall) $ min (upperBound fbc :: CN MPBall) cutoff+ leftMonotone = V.foldl' (||) False $ V.map (!/=! 0) dfa'+ rightMonotone = V.foldl' (||) False $ V.map (!/=! 0) dfb'+ boxes = + case (leftMonotone || fa' !>! cutoff', rightMonotone || fb' !>! cutoff') of+ (True, True) -> []+ (True, False) -> [SearchBox b fb']+ (False, True) -> [SearchBox a fa']+ (False, False) -> [SearchBox a fa', SearchBox b fb']+ in+ (cutoff', boxes)++-- Precondition: f and g must have the same domain+maxBoxFunGreaterThanN :: BoxFun -> BoxFun -> CN Rational -> Precision -> Bool+maxBoxFunGreaterThanN f g n initialPrecision =+ case Box.getEndpoints fbox == Box.getEndpoints gbox of+ CertainTrue ->+ checkMaxAboveN f g ||+ (Box.width fboxp !>! cutoff && Box.width gboxp !>! cutoff) &&+ let+ newBoxes = Box.fullBisect fboxp++ updateDomain z = BoxFun (dimension z) (bf_eval z)++ checkBoxes [] = True+ checkBoxes (box : boxes) = + if checkMaxAboveN (updateDomain f box) (updateDomain g box) + then checkBoxes boxes+ else maxBoxFunGreaterThanN f' g' n initialPrecision && checkBoxes boxes+ where+ f' = updateDomain f box+ g' = updateDomain g box++ in+ checkBoxes newBoxes++ _ ->+ trace "Domain of f not equal to domain of g"+ False+ where+ cutoff = 1/2^10++ fbox = domain f+ fboxp = setPrecision initialPrecision fbox ++ gbox = domain g+ gboxp = setPrecision initialPrecision gbox++ checkMaxAboveN h i = applyMinimum h !>! n || applyMinimum i !>! n++
+ src/AERN2/BoxFun/TestFunctions.hs view
@@ -0,0 +1,425 @@+module AERN2.BoxFun.TestFunctions where++import MixedTypesNumPrelude+import AERN2.MP.Ball+import qualified Data.List as List++import AERN2.AD.Differential+import AERN2.BoxFun.Type+import AERN2.Linear.Vector.Type ((!), Vector)+import qualified AERN2.Linear.Vector.Type as V++fromListDomain :: [(Rational, Rational)] -> Vector (CN MPBall)+fromListDomain [] = V.empty+fromListDomain [x] = V.singleton $ fromEndpointsAsIntervals (cn $ mpBallP (prec 53) $ (fst x)) (cn $ mpBallP (prec 53) (snd x))+fromListDomain (x:xs) = V.cons (fromEndpointsAsIntervals (cn $ mpBallP (prec 53) $ (fst x)) (cn $ mpBallP (prec 53) (snd x))) (fromListDomain xs)++symmetricDomain :: Integer -> Rational -> Rational -> Vector (CN MPBall)+symmetricDomain n l r = + V.map (\_ -> fromEndpointsAsIntervals (cn $ mpBallP (prec 53) $ l) (cn $ mpBallP (prec 53) r)) $ V.enumFromTo 1 n++griewank :: Integer -> BoxFun+griewank n =+ BoxFun+ n+ (\v ->+ let+ ord = order (v ! 0)+ sm = List.foldl' (+) (differential ord (cn $ mpBall 0)) [let x = (v ! k) in x^2 | k <- [0 .. V.length v - 1]]+ prd = List.foldl' (*) (differential ord (cn $ mpBall 1)) [cos $ (v ! k) / (sqrt $ 1 + mpBallP p k) | k <- [0 .. V.length v - 1]]+ p = getPrecision v+ in+ 1 + ((setPrecision p $ mpBall 1)/(setPrecision p $ mpBall 4000)) * sm - prd+ )+ (symmetricDomain n (-600.0) 600.0)++rosenbrock :: BoxFun+rosenbrock =+ BoxFun+ 2+ (\v ->+ let+ p = getPrecision v+ x = v ! 0+ y = v ! 1+ a = mpBallP p 1+ b = mpBallP p 100+ amx = a - x+ ymxs = y - x*x+ in+ amx^2 + b*ymxs^2+ )+ (symmetricDomain 2 (-1.2) 1.2)++himmelblau :: BoxFun+himmelblau =+ BoxFun+ 2+ (\v ->+ let+ x = v ! 0+ y = v ! 1+ a = (x^2 + y - 11)+ b = (x + y^2 - 7)+ in+ a^2 + b^2+ )+ (symmetricDomain 2 (-600.0) 600.0)++shekel :: BoxFun+shekel = + BoxFun+ 2+ (\v ->+ let+ x = v ! 0+ y = v ! 1+ c0 = 1; a00 = 43; a01 = 23+ c1 = 17; a10 = 6; a11 = 9+ in+ - 1/(c0 + (x - a00)^2 + (y - a01)^2)+ - 1/(c1 + (x - a10)^2 + (y - a11)^2)+ )+ (symmetricDomain 2 (-600.0) 600.0)++siam4 :: BoxFun+siam4 = + BoxFun+ 2+ (\v ->+ let+ x = v ! 0+ y = v ! 1+ in+ exp(sin(50 * x)) + sin(60 * exp y) + + sin(70 * sin(x)) + sin(sin(80 * y)) + - sin(10 * (x + y)) + + (x^2 + y^2) / 4+ )+ (symmetricDomain 2 (-10.0) 10.0)++ratz4 :: BoxFun+ratz4 =+ BoxFun+ 2+ (\v ->+ let + x = v ! 0+ y = v ! 1+ xs = x^2+ ys = y^2+ in + sin(xs + 2 * ys) * exp (-xs - ys)+ )+ (symmetricDomain 2 (-3.0) 3.0)++-- global minimum at bukin(-10, 1) ~ 0+-- bukin :: BoxFun+-- bukin =+-- BoxFun+-- 2+-- (\v ->+-- let+-- x = v!0+-- y = v!1+-- in+-- 100 * sqrt (abs (y - x^2 / 100)) + abs(x + 10) / 100+-- )+-- (fromListDomain [(-15.0, 5.0), (-3.0, 3.0)])++ackley :: BoxFun+ackley =+ BoxFun+ 2+ (\v ->+ let+ x = v!0+ y = v!1+ p = getPrecision v+ pi = piBallP p+ in+ -20 * exp(sqrt((x^2 + y^2) / 2) / (-5)) - exp((cos (2 * pi * x) + cos(2 * pi * y)) / 2) + exp(mpBallP p 1) + 20+ )+ (symmetricDomain 2 (-5.0) 5.0)++-- eggholder :: BoxFun+-- eggholder =+-- BoxFun+-- 2+-- (\v ->+-- let+-- x = v!0+-- y = v!1+-- in+-- -(y + 47) * sin (sqrt (abs (x / 2 + (y + 47)))) - x * sin (sqrt (abs (x - (y + 47))))+-- )+-- (symmetricDomain 2 (-512.0) 512.0)+ +-- heron :: BoxFun+-- heron = +-- BoxFun+-- 2+-- (\v ->+-- let+-- x = v!0+-- y = v!1+-- p = getPrecision v+-- eps = 1/2^(23)+-- i = 2+-- in+-- max+-- min+-- max+-- (y - sqrt x)+-- ((sqrt x - y) - (mpBallP p 1/2)^(2^(i-1)) - 6 * eps * (i-1))+-- max+-- (sqrt x - y)+-- (- (sqrt x - y) - (mpBallP p 1/2)^(2^(i-1)) - 6 * eps * (i-1))+-- min+-- max+-- ((y + x/y)/2 - sqrt x)+-- (- (sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + 6 * eps * (i-1))+-- max+-- (sqrt x - (y+x/y)/2)+-- ((sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + 6 * eps * (i-1))++-- )+-- (fromListDomain [(0.5, 2.0), (0.8, 1.8)])++-- max (min (max 1p 1q) (max 2p 2q)) (min (max 3p 3q) (max 4p 4q))++i :: Integer+i = 3++heron1p :: BoxFun+heron1p = + BoxFun+ 2+ (\v ->+ let+ x = v!0+ y = v!1+ _p = getPrecision v+ _eps = 1/2^(23)+ in+ (y - sqrt x)+ )+ (fromListDomain [(0.5, 2.0), (0.8, 1.8)])++ +heron1q :: BoxFun+heron1q = + BoxFun+ 2+ (\v ->+ let+ x = v!0+ y = v!1+ p = getPrecision v+ eps = 1/2^(23)+ in+ ((sqrt x - y) - (mpBallP p 1/2)^(2^(i-1)) - (mpBallP p 6) * eps * (i-1))++ )+ (fromListDomain [(0.5, 2.0), (0.8, 1.8)])++ +heron2p :: BoxFun+heron2p = + BoxFun+ 2+ (\v ->+ let+ x = v!0+ y = v!1+ _p = getPrecision v+ _eps = 1/2^(23)+ in+ (sqrt x - y)+ )+ (fromListDomain [(0.5, 2.0), (0.8, 1.8)])++ +heron2q :: BoxFun+heron2q = + BoxFun+ 2+ (\v ->+ let+ x = v!0+ y = v!1+ p = getPrecision v+ eps = 1/2^(23)+ in+ (- (sqrt x - y) - (mpBallP p 1/2)^(2^(i-1)) - (mpBallP p 6) * eps * (i-1))++ )+ (fromListDomain [(0.5, 2.0), (0.8, 1.8)])++ +heron3p :: BoxFun+heron3p = + BoxFun+ 2+ (\v ->+ let+ x = v!0+ y = v!1+ _p = getPrecision v+ _eps = 1/2^(23)+ in+ ((y + x/y)/2 - sqrt x)+ )+ (fromListDomain [(0.5, 2.0), (0.8, 1.8)])+ +heron3q :: BoxFun+heron3q = + BoxFun+ 2+ (\v ->+ let+ x = v!0+ y = v!1+ p = getPrecision v+ eps = 1/2^(23)+ in+ (- (sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + (mpBallP p 6) * eps * (i-1))+ )+ (fromListDomain [(0.5, 2.0), (0.8, 1.8)])++ +heron3p2 :: BoxFun+heron3p2 = + BoxFun+ 2+ (\v ->+ let+ x = v!0+ y = v!1+ _p = getPrecision v+ _eps = 1/2^(23)+ _i = 4+ in+ ((y + x/y)/2 - sqrt x)+ )+ (fromListDomain [(0.633758544921875, 0.63385009765625), (0.79999999999999982236431605997495353221893310546875, 0.80006103515624982239142111428709114306911942549049854278564453125)])+ +heron3q2 :: BoxFun+heron3q2 = + BoxFun+ 2+ (\v ->+ let+ x = v!0+ y = v!1+ p = getPrecision v+ eps = 1/2^(23)+ _i = 4+ in+ (- (sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + (mpBallP p 6) * eps * (i-1))+ )+ (fromListDomain [(0.633758544921875, 0.63385009765625), (0.79999999999999982236431605997495353221893310546875, 0.80006103515624982239142111428709114306911942549049854278564453125)])++heron4p :: BoxFun+heron4p = + BoxFun+ 2+ (\v ->+ let+ x = v!0+ y = v!1+ _p = getPrecision v+ _eps = 1/2^(23)+ in+ (sqrt x - (y+x/y)/2)+ )+ (fromListDomain [(0.5, 2.0), (0.8, 1.8)])++ +heron4q :: BoxFun+heron4q = + BoxFun+ 2+ (\v ->+ let+ x = v!0+ y = v!1+ p = getPrecision v+ eps = 1/2^(23)+ in+ ((sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + (mpBallP p 6) * eps * (i-1))++ )+ (fromListDomain [(0.5, 2.0), (0.8, 1.8)])++-- heronFull :: BoxFun+-- heronFull =+-- BoxFun+-- 3+-- (\v ->+-- let+-- x = v!0+-- y = v!1+-- i = v!2+-- p = getPrecision v+-- eps = 1/2^(23)+-- in+-- ((sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + (mpBallP p 6) * eps * (i-1))++-- )+-- (fromListDomain [(0.5, 2.0), (0.8, 1.8), (1.0, 5.0)])+++mxp1 :: BoxFun+mxp1 =+ BoxFun+ 2+ (\v ->+ let+ x = v!0+ -- y = v!1+ in+ -x+1++ )+ (fromListDomain [(0.0, 2.0), (0.0, 2.0)])++xm1 :: BoxFun+xm1 =+ BoxFun+ 2+ (\v ->+ let+ x = v!0+ -- y = v!1+ in+ x-1++ )+ (fromListDomain [(0.0, 2.0), (0.0, 2.0)])++xe2p1 :: BoxFun+xe2p1 =+ BoxFun+ 2+ (\v ->+ let+ x = v!0+ -- y = v!1+ in+ x^2+1+ )+ (fromListDomain [(0.0, 2.0), (0.0, 2.0)])++xe2m1 :: BoxFun+xe2m1 =+ BoxFun+ 2+ (\v ->+ let+ x = v!0+ -- y = v!1+ in+ x^2-1+ )+ (fromListDomain [(0.0, 2.0), (0.0, 2.0)])
+ src/AERN2/BoxFun/Type.hs view
@@ -0,0 +1,155 @@+module AERN2.BoxFun.Type where++import MixedTypesNumPrelude++import AERN2.Linear.Vector.Type (Vector, (!))+import qualified AERN2.Linear.Vector.Type as V+import AERN2.MP.Ball+import AERN2.AD.Differential+import qualified AERN2.Linear.Matrix.Type as M+import AERN2.Util.Util+import AERN2.BoxFun.Box+import Numeric.CollectErrors++data BoxFun =+ BoxFun+ {+ dimension :: Integer+ , bf_eval :: Vector (Differential (CN MPBall)) -> Differential (CN MPBall)+ , domain :: Vector (CN MPBall)+ }++instance Show BoxFun where+ show (BoxFun d _f b) = show d ++ " dimensional BoxFun, domain: " ++ show b++boundaryRestrictions :: BoxFun -> [BoxFun]+boundaryRestrictions (BoxFun d ev dom) =+ concat+ [+ [+ BoxFun + (d - 1)+ (+ \v ->+ ev $ V.map (\j -> if j == i then (setPrecision $ getPrecision v) $ differential 2 $ upperBound (dom ! i) else if j < i then v ! j else v ! (j - 1)) $ V.enumFromTo 0 (d - 1)+ ) + (V.map (\j -> if j >= i then dom ! (j + 1) else dom ! j) $ V.enumFromTo 0 (d - 2))+ ,+ BoxFun + (d - 1)+ (+ \v ->+ ev $ V.map (\j -> if j == i then (setPrecision $ getPrecision v) $ differential 2 $ lowerBound (dom ! i) else if j < i then v ! j else v ! (j - 1)) $ V.enumFromTo 0 (d - 1)+ )+ (V.map (\j -> if j >= i then dom ! (j + 1) else dom ! j) $ V.enumFromTo 0 (d - 2))+ ]+ |+ i <- [0 .. d - 1]+ ]++valueGradientHessian :: BoxFun -> Vector (CN MPBall) -> (CN MPBall, Vector (CN MPBall), M.Matrix (CN MPBall))+valueGradientHessian (BoxFun d e _) v =+ (value, grad, hess)+ where+ vgh i j = e (w i j)++ triangle = + V.map (\i-> V.map (\j -> vgh i j) $ V.enumFromTo 0 i) $ V.enumFromTo 0 (d - 1)++ value = diff_x $ (triangle ! 0) ! 0+ grad = V.map (\i -> diff_dx $ (triangle ! i) ! 0) $ V.enumFromTo 0 (d - 1)+ hess = M.create d d (\i j -> diff_d2x $ if i > j then (triangle ! i) ! j else (triangle ! j) ! i)++ w i j = V.imap (\k x -> OrderTwo x (delta i k) (delta j k) (pure $ mpBall 0)) v+ delta :: Integer -> Integer -> CN MPBall+ delta i k = if i == k then (cn $ mpBall 1) else (cn $ mpBall 0)++valueGradient :: BoxFun -> Vector (CN MPBall) -> (CN MPBall, Vector (CN MPBall))+valueGradient (BoxFun d e _) v =+ aux (d - 1) [] (pure $ mpBall 0)+ where+ tangent k = + V.imap (\i x -> OrderOne x (delta i k)) v+ valgrad k =+ let+ etk = e (tangent k)+ in + (diff_x etk, diff_dx etk)+ aux k ret val =+ if k < 0 then+ (val, V.fromList ret)+ else + let+ (val2,g) = valgrad k+ in+ aux (k - 1) (g : ret) val2+ delta :: Integer -> Integer -> CN MPBall+ delta i k = if i == k then (cn $ mpBall 1) else (cn $ mpBall 0)++apply :: BoxFun -> Vector (CN MPBall) -> CN MPBall+apply (BoxFun _d e _) v = + diff_x (e v')+ where+ v' = V.map (\x -> differential 0 x) v++applyMinimum :: BoxFun -> CN MPBall+applyMinimum h = fst $ endpointsAsIntervals (apply h (domain h))++applyMinimumOnBox :: BoxFun -> Vector (CN MPBall) -> CN MPBall+applyMinimumOnBox h hbox = fst $ endpointsAsIntervals (apply h hbox)++applyMaximum :: BoxFun -> CN MPBall+applyMaximum h = snd $ endpointsAsIntervals (apply h (domain h))++applyMaximumOnBox :: BoxFun -> Vector (CN MPBall) -> CN MPBall+applyMaximumOnBox h hbox = snd $ endpointsAsIntervals (apply h hbox)++gradient :: BoxFun -> Vector (CN MPBall) -> Vector (CN MPBall)+gradient (BoxFun d e _) v =+ aux (d - 1) []+ where+ tangent k = + V.imap (\i x -> OrderOne x (delta i k)) v+ grad k =+ -- diff_dx $ e (tangent k)+ case e (tangent k) of+ OrderZero _ -> noValueNumErrorCertain (NumError "details")+ val -> diff_dx $ val+ aux k ret =+ if k < 0 then+ V.fromList ret+ else + aux (k - 1) (grad k : ret)+ delta :: Integer -> Integer -> CN MPBall+ delta i k = if i == k then (cn $ mpBall 1) else (cn $ mpBall 0)++hessian :: BoxFun -> Vector (CN MPBall) -> M.Matrix (CN MPBall)+hessian (BoxFun d e _) v = + M.create d d a+ where+ a i j = diff_d2x $ e (w i j)+ w i j = V.imap (\k x -> OrderTwo x (delta i k) (delta j k) (pure $ mpBall 0)) v+ delta :: Integer -> Integer -> CN MPBall+ delta i k = if i == k then (cn $ mpBall 1) else (cn $ mpBall 0)++jacobian :: [BoxFun] -> Vector (CN MPBall) -> M.Matrix (CN MPBall)+jacobian fs v =+ M.create (length fs) highestDiminseionInFs a+ where+ highestDiminseionInFs = maximum (map dimension fs)++ a i j = diff_dx $ bf_eval (fs!!i) (w j)+ w j = V.imap (\k x -> OrderOne x (delta j k)) v+ delta :: Integer -> Integer -> CN MPBall+ delta i k = if i == k then (cn $ mpBall 1) else (cn $ mpBall 0)++gradientUsingGradient :: BoxFun -> Box -> Box+gradientUsingGradient f v =+ V.zipWith fromEndpointsAsIntervals lowerBounds upperBounds+ where+ lowerBounds = firstDerivatives - secondDerivatives * V.map (fmap (mpBall . radius)) v+ upperBounds = firstDerivatives + secondDerivatives * V.map (fmap (mpBall . radius)) v+ firstDerivatives = head $ M.rows $ jacobian [f] (centre v)+ secondDerivatives = fmap abs $ hessian f v+ p = getPrecision v+
+ src/AERN2/Linear/Matrix/Inverse.hs view
@@ -0,0 +1,126 @@+module AERN2.Linear.Matrix.Inverse +(inverse) +where++import MixedTypesNumPrelude+import qualified Numeric.CollectErrors as CN++import AERN2.Linear.Matrix.Type+import qualified AERN2.Linear.Vector.Type as V+import Data.Maybe+import AERN2.MP.Float+import AERN2.MP.Dyadic+import AERN2.MP++inverse :: Matrix (CN MPBall) -> Maybe (Matrix (CN MPBall))+inverse m = + if CN.hasError m+ || (isNothing maybeY)+ || (not $ nerr !<! 1) then+ Nothing+ else+ Just $ aux x0+ where+ cm = Matrix (width m) (V.map (\x -> mpFloat $ centre $ unCN x) $ entries m)+ idM = identity (width m) (width m) :: Matrix (CN MPBall)+ maybeY = inverseGauss cm+ y = fmap (cn . mpBall . dyadic) $ fromJust maybeY+ yid = y * idM+ nye = V.map V.inftyNorm (columns yid)+ err = idM - y * m+ nerr = inftyNorm err+ ui = cn $ fromEndpointsAsIntervals (mpBall $ -1) (mpBall 1)+ x0 = create (width m) (width m) (\_ j -> (ui * (nye V.! j)) / (1 - nerr))+ aux x = let x' = it x in if getAccuracy x' <= getAccuracy x then x' else aux x'+ it x = intersect x $ yid + err * x++inverseGauss :: Matrix MPFloat -> Maybe (Matrix MPFloat)+inverseGauss m = + aux 0 m inv0+ where+ inv0 = identity (width m) (width m) :: Matrix MPFloat+ aux j n inv = + if j == width m then+ Just inv+ else+ let+ n' = pivot n j n+ inv' = pivot n j inv+ a = get n' j j + ia = ((mpFloat 1)/.a)+ n'' = multiplyRow n' j ia + inv'' = multiplyRow inv' j ia + elim i l = + if i /= j then+ let+ b = get l i j+ in+ addRows l (-b) j i+ else + l+ n''' = + foldr elim n'' [0 .. (height m) - 1]+ inv''' = + foldr elim inv'' [0 .. (height m) - 1]+ in + if a /= mpFloat 0 then+ aux (j + 1) n''' inv'''+ else + Nothing+++largestRowIndex :: Matrix MPFloat -> Integer -> Integer+largestRowIndex m j =+ integer $ aux (j + 1) j (get m j j)+ where+ h = height m+ aux k i x =+ if k == h then+ i+ else + let+ x' = get m i j+ in + if x' > x then+ aux (k + 1) k x'+ else + aux (k + 1) i x+ +pivot :: Matrix MPFloat -> Integer -> Matrix MPFloat -> Matrix MPFloat+pivot m j n = + swapRows n j i+ where+ i = largestRowIndex m j++multiplyRow :: Matrix MPFloat -> Integer -> MPFloat -> Matrix MPFloat+multiplyRow m i a = + imap h m+ where+ h k _ x = + if k == i then + a *. x+ else + x++{- add a*(row i0) to row i1 -}+addRows :: Matrix MPFloat -> MPFloat -> Integer -> Integer -> Matrix MPFloat+addRows m a i0 i1 = + imap h m+ where+ h k l x = + if k == i1 then + x +. a *. get m i0 l+ else + x++swapRows :: Matrix MPFloat -> Integer -> Integer -> Matrix MPFloat+swapRows m i0 i1 =+ imap h m+ where+ h k l x = + if k == i0 then + get m i1 l+ else if k == i1 then+ get m i0 l+ else + x
+ src/AERN2/Linear/Matrix/Type.hs view
@@ -0,0 +1,132 @@+module AERN2.Linear.Matrix.Type where++import qualified Prelude as P+import MixedTypesNumPrelude+import qualified Numeric.CollectErrors as CN+import AERN2.Linear.Vector.Type (Vector, (!))+import qualified AERN2.Linear.Vector.Type as V+import qualified Data.List as List++import AERN2.MP.Ball++data (Matrix a) = + Matrix+ {+ width :: Integer,+ entries :: Vector a+ } deriving (Show)++height :: Matrix a -> Integer+height (Matrix w e) = + (V.length e) `P.div` w++get :: Matrix a -> Integer -> Integer -> a+get m i j =+ entries m ! (i * (width m) + j)++identity :: (HasIntegers a) => Integer -> Integer -> Matrix a+identity m n = + diag m n (convertExactly 1)++diag :: (HasIntegers a) => Integer -> Integer -> a -> Matrix a+diag m n x = + create m n (\i j -> if i == j then x else (convertExactly 0))++rows :: Matrix a -> [Vector a]+rows m@(Matrix w e) = + [V.slice (i*w) w e| i <- [0 .. height m - 1]]++columns :: Matrix a -> Vector (Vector a)+columns m = + V.map (\j -> V.map (\i -> get m i j) $ V.enumFromTo 0 (height m - 1)) $ V.enumFromTo 0 (width m - 1)++create :: Integer -> Integer -> (Integer -> Integer -> a) -> Matrix a+create m n f =+ Matrix n $ V.map (\x -> f (i x) (j x)) $ V.enumFromTo 0 (m*n - 1)+ where+ j x = x `mod` n+ i x = (x - j x) `P.div` n++imap :: (Integer -> Integer -> a -> a) -> Matrix a -> Matrix a+imap f (Matrix w ents) =+ Matrix w (V.imap g ents)+ where+ j x = x `mod` w+ i x = (x - j x) `P.div` w+ g k x = f (i k) (j k) x++instance CanIntersectAsymmetric (Matrix (CN MPBall)) (Matrix (CN MPBall)) where+ type IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall)) = Matrix (CN MPBall)+ intersect (Matrix w0 v0) (Matrix _w1 v1) =+ Matrix w0 $ V.zipWith intersect v0 v1++inftyNorm :: (CanAddSameType a, CanSubSameType a, CanAbsSameType a, HasIntegers a, CanMinMaxSameType a) => Matrix a -> a+inftyNorm (m :: Matrix a) = + -- TODO: could be optimised.+ List.foldl' max (convertExactly 0 :: a)+ [+ V.foldl' (+) (convertExactly 0 :: a) $ V.map abs r + |+ r <- rows m+ ]++instance Functor Matrix where+ fmap h m = + Matrix (width m) (V.map h (entries m))++instance + (CanAddSameType a, CanMulSameType a, HasIntegers a) =>+ CanMulAsymmetric (Matrix a) (Matrix a)+ where+ type MulType (Matrix a) (Matrix a) = Matrix a+ mul m0 m1 = + create (height m0) (width m1) (aux 0 (convertExactly 0))+ where+ aux k sm i j = + if k == width m0 then + sm+ else + aux (k + 1) (sm + (get m0 i k) * (get m1 k j)) i j+++instance + (CanAddSameType a) =>+ CanAddAsymmetric (Matrix a) (Matrix a)+ where+ type AddType (Matrix a) (Matrix a) = Matrix a+ add (Matrix w e) (Matrix _ e') =+ Matrix w (e + e')++instance + (CanSubSameType a) =>+ CanSub (Matrix a) (Matrix a)+ where+ type SubType (Matrix a) (Matrix a) = Matrix a+ sub (Matrix w e) (Matrix _ e') =+ Matrix w (e - e')+ ++instance + (CanAddSameType a, CanMulSameType a, HasIntegers a) =>+ CanMulAsymmetric (Matrix a) (Vector a)+ where+ type MulType (Matrix a) (Vector a) = Vector a+ mul m@(Matrix _w _e) v =+ V.fromList [r * v| r <- rows m]++instance + (HasAccuracy a, HasPrecision a) => HasAccuracy (Matrix a)+ where+ getAccuracy m = + V.foldl' max NoInformation $ V.map getAccuracy (entries m)++instance + (HasPrecision a) => HasPrecision (Matrix a)+ where+ getPrecision m = + V.foldl' max (prec 2) $ V.map getPrecision (entries m)++instance + (CN.CanTestErrorsPresent a) => CN.CanTestErrorsPresent (Matrix a)+ where+ hasError m = V.foldl' (||) False $ V.map (CN.hasError) (entries m)
+ src/AERN2/Linear/Vector/Type.hs view
@@ -0,0 +1,176 @@+module AERN2.Linear.Vector.Type where++import Control.Monad.ST+import Data.STRef+import MixedTypesNumPrelude hiding (length)+import qualified Data.Vector as V+import qualified Data.Vector.Generic.Mutable as M+import AERN2.MP.Precision+import AERN2.MP.Ball+import qualified Prelude as P++type (Vector a) = V.Vector a++(+++) :: Vector a -> Vector a -> Vector a+(+++) = (V.++)++drop :: Int -> Vector a -> Vector a +drop = V.drop++take :: Int -> Vector a -> Vector a +take = V.take++empty :: Vector a+empty = V.empty++singleton :: a -> Vector a+singleton = V.singleton++cons :: a -> Vector a -> Vector a+cons = V.cons++fromList :: [a] -> Vector a+fromList = V.fromList++map :: (a -> b) -> Vector a -> Vector b+map = V.map++imap :: (Integer -> a -> b) -> Vector a -> Vector b+imap h = V.imap (\i x -> h (integer i) x)++enumFromTo :: Enum a => a -> a -> Vector a+enumFromTo = V.enumFromTo++slice :: Integer -> Integer -> Vector a -> Vector a+slice i j = V.slice (int i) (int j)++foldl' :: (b -> a -> b) -> b -> Vector a -> b+foldl' = V.foldl'++zipWith :: (a -> b -> c) -> Vector a -> Vector b -> Vector c+zipWith = V.zipWith++(!) :: Vector a -> Integer -> a+(!) v i = (V.!) v (int i)++length :: Vector a -> Integer+length = integer . V.length++intLength :: Vector a -> Int +intLength = V.length++inftyNorm :: (HasIntegers a, CanMinMaxSameType a) => Vector a -> a+inftyNorm (v :: Vector a) =+ V.foldl' max (convertExactly 0 :: a) v++find :: (a -> Bool) -> Vector a -> Maybe a +find = V.find++elem :: P.Eq a => a -> Vector a -> Bool+elem = V.elem++toList :: Vector a -> [a] +toList = V.toList++zip :: Vector a -> Vector b -> Vector (a, b) +zip = V.zip++null :: Vector a -> Bool+null = V.null++any :: (a -> Bool) -> V.Vector a -> Bool+any = V.any++all :: (a -> Bool) -> V.Vector a -> Bool+all = V.all+instance + (HasAccuracy a, HasPrecision a) => HasAccuracy (Vector a)+ where+ getAccuracy v = + V.foldl' max NoInformation $ V.map getAccuracy v++instance + (HasPrecision a) => HasPrecision (Vector a)+ where+ getPrecision v = + if V.null v then + (prec 2) + else + getPrecision $ v ! 0 -- TODO: safe? Alternative: V.foldl' max (prec 2) $ V.map getPrecision v ++instance + (CanSetPrecision a) => CanSetPrecision (Vector a)+ where+ setPrecision p = V.map (setPrecision p)++instance + (CanAddSameType a) =>+ CanAddAsymmetric (Vector a) (Vector a)+ where+ type AddType (Vector a) (Vector a) = Vector a+ add v w =+ runST $+ do+ mv <- M.new (intLength v)+ aux mv 0+ V.freeze mv+ where+ lth = length v+ aux :: (V.MVector s a) -> Integer -> (ST s ())+ aux mv k = + if k == lth then+ return ()+ else+ do + M.write mv (int k) (v ! k + w ! k)+ aux mv (k + 1)++instance + (CanSubSameType a) =>+ CanSub (Vector a) (Vector a)+ where+ type SubType (Vector a) (Vector a) = Vector a+ sub v w =+ runST $+ do+ mv <- M.new (intLength v)+ aux mv 0+ V.freeze mv+ where+ lth = length v+ aux :: (V.MVector s a) -> Integer -> (ST s ())+ aux mv k = + if k == lth then+ return ()+ else+ do + M.write mv (int k) (v ! k - w ! k)+ aux mv (k + 1)+++instance + (CanAddSameType a, CanMulSameType a, HasIntegers a) =>+ CanMulAsymmetric (Vector a) (Vector a) + where+ type MulType (Vector a) (Vector a) = a+ mul v w =+ runST $+ do+ sum <- newSTRef (convertExactly 0)+ aux sum 0+ readSTRef sum+ where+ lth = length v+ aux :: (STRef s a) -> Integer -> (ST s ())+ aux sum k = + if k == lth then+ return ()+ else+ do + modifySTRef sum (\x -> x + (v ! k) * (w ! k))+ aux sum (k + 1)++instance + CanMulAsymmetric (CN MPBall) (Vector (CN MPBall)) where+ type MulType (CN MPBall) (Vector (CN MPBall)) = Vector (CN MPBall)+ mul x v = V.map (\y -> x * y) v
+ src/AERN2/Util/Util.hs view
@@ -0,0 +1,14 @@+module AERN2.Util.Util where++import MixedTypesNumPrelude++import AERN2.MP.Ball++upperBound :: (IsInterval i) => i -> i+upperBound = endpointRAsInterval++lowerBound :: (IsInterval i) => i -> i+lowerBound = endpointLAsInterval++intersectCN :: CN MPBall -> CN MPBall -> CN MPBall+intersectCN = intersectCNMPBall