AERN-RnToRm (empty) → 0.3.0
raw patch · 26 files changed
+6400/−0 lines, 26 filesdep +AERN-Realdep +basedep +binarysetup-changed
Dependencies added: AERN-Real, base, binary, containers
Files
- AERN-RnToRm.cabal +123/−0
- ChangeLog +4/−0
- LICENCE +30/−0
- Setup.lhs +3/−0
- src/Data/Number/ER/RnToRm.hs +100/−0
- src/Data/Number/ER/RnToRm/Approx.hs +303/−0
- src/Data/Number/ER/RnToRm/Approx/DomEdges.hs +487/−0
- src/Data/Number/ER/RnToRm/Approx/DomTransl.hs +495/−0
- src/Data/Number/ER/RnToRm/Approx/PieceWise.hs +522/−0
- src/Data/Number/ER/RnToRm/Approx/Tuple.hs +272/−0
- src/Data/Number/ER/RnToRm/BisectionTree.hs +665/−0
- src/Data/Number/ER/RnToRm/BisectionTree/Integration.hs +278/−0
- src/Data/Number/ER/RnToRm/BisectionTree/Path.hs +134/−0
- src/Data/Number/ER/RnToRm/DefaultRepr.hs +64/−0
- src/Data/Number/ER/RnToRm/TestingDefs.hs +72/−0
- src/Data/Number/ER/RnToRm/UnitDom/Approx.hs +92/−0
- src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs +589/−0
- src/Data/Number/ER/RnToRm/UnitDom/Base.hs +348/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs +93/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs +279/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs +293/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs +455/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs +190/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs +226/−0
- src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs +169/−0
- tests/Demo.hs +114/−0
+ AERN-RnToRm.cabal view
@@ -0,0 +1,123 @@+Name: AERN-RnToRm+Version: 0.3.0+Cabal-Version: >= 1.2+Build-Type: Simple+License: BSD3+License-File: LICENCE+Author: Michal Konecny (Aston University)+Copyright: (c) 2007-2008 Michal Konecny+Maintainer: mik@konecny.aow.cz+Stability: experimental+Category: Data, Math+Synopsis: polynomial function enclosures (PFEs) approximating exact real functions+Tested-with: GHC ==6.8.2+Description:+ AERN-RnToRm provides+ datatypes and abstractions for approximating functions+ of type @D -> R^m@ where @D@ is a bounded interval in @R^n@+ with non-empty interior.+ .+ Abstractions are provided via 4 type classes:+ .+ * ERUnitFnBase: + generalises polynomials with floating point coefficients.+ (/Not exported here, used only internally./)+ .+ * ERFnApprox:+ generalises functions enclosures on a certain unspecified domain.+ .+ * ERUnitFnApprox (extends ERFnApprox): generalises function graph enclosures+ on the domain @[-1,1]^n@.+ (/Not exported here, used only internally./)+ .+ * ERFnDomApprox (extends ERFnApprox):+ generalises function enclosures over a specified and queriable domain box+ (instance of class DomainBox).+ .+ At all levels, all field operations are supported as well as+ some elementary operations, namely exp, sin and cos.+ Log and sqrt are planned to be added soon. + .+ Implementations of ERUnitFnBase:+ .+ * ERChebPoly + .+ By using the Chebyshev basis on domain @[-1,1]^n@, + we gain simple and optimally rounding degree reduction + as well as relatively simple handling of rounding + in other operations.+ .+ Implementations of ERUnitFnApprox:+ . + * ERFnInterval+ .+ Implementations of ERFnDomApprox:+ .+ * ERFnDomTranslApprox: + builds a basic implementation + using an instance of ERUnitFnApprox.+ .+ * ERFnTuple: + extends another implementation of ERFnDomApprox + to work with lists of functions simultaneously.+ .+ * ERFnDomEdgesApprox: + separately enclose a function on its domain box + as well as on all the domain's hyper-edges + (including the corners) using+ another implementation of ERFnDomApprox.+ .+ * ERFnPiecewise: + allows the domain box to be bisected + to an arbitrary finite depth + and uses another implementation of ERFnDomApprox + to approximate the function on each segment. + .+ Simple examples of usage can be found in tests/Demo.hs.+ +Extra-source-files:+ ChangeLog tests/Demo.hs++Flag containers-in-base+ Default: False++Library+ hs-source-dirs: src+ if flag(containers-in-base)+ Build-Depends:+ base < 3, binary >= 0.4, AERN-Real == 0.9.6+ else+ Build-Depends:+ base >= 3, containers, binary >= 0.4, AERN-Real == 0.9.6+ Exposed-modules:+ Data.Number.ER.RnToRm,+ Data.Number.ER.RnToRm.BisectionTree.Path,+ Data.Number.ER.RnToRm.BisectionTree.Integration,+ Data.Number.ER.RnToRm.BisectionTree,+ Data.Number.ER.RnToRm.DefaultRepr,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Integration,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds,+ Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom,+ Data.Number.ER.RnToRm.UnitDom.Base,+ Data.Number.ER.RnToRm.UnitDom.Approx.Interval,+ Data.Number.ER.RnToRm.UnitDom.Approx,+ Data.Number.ER.RnToRm.Approx.DomTransl,+ Data.Number.ER.RnToRm.Approx.PieceWise,+ Data.Number.ER.RnToRm.Approx.DomEdges,+ Data.Number.ER.RnToRm.Approx.Tuple,+ Data.Number.ER.RnToRm.Approx,+ Data.Number.ER.RnToRm.TestingDefs + + Extensions: + CPP,+ DeriveDataTypeable,+ FlexibleContexts,+ FlexibleInstances,+ FunctionalDependencies,+ MultiParamTypeClasses,+ UndecidableInstances+
+ ChangeLog view
@@ -0,0 +1,4 @@+0.3.0: 7 August 2008+ * initial release of AERN-RnToRm after one year of work and two successful + internal applications+
+ LICENCE view
@@ -0,0 +1,30 @@+Copyright (c) 2007-2008 Michal Konecny, Amin Farjudian, Jan Duracz++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.++3. Neither the name of the author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE CONTRIBUTORS ``AS IS'' AND ANY EXPRESS+OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,+STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE+POSSIBILITY OF SUCH DAMAGE.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ src/Data/Number/ER/RnToRm.hs view
@@ -0,0 +1,100 @@+{-|+ Module : Data.Number.ER.RnToRm+ Description : overview of AERN-RnToRm+ Copyright : (c) Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : non-portable (requires fenv.h)++ This module bundles some of the most important functionality+ of the AERN-RnToRm package. + It is intended to be imported *qualified*.++ AERN-RnToRm provides+ datatypes and abstractions for approximating functions+ of type @D -> R^m@ where @D@ is a bounded interval in @R^n@+ with non-empty interior.+ + Abstractions are provided via 4 type classes:+ + * 'UFB.ERUnitFnBase': + generalises polynomials with floating point coefficients.+ (/Not exported here, used only internally./)+ + * 'ERFnApprox': + generalises functions enclosures on a certain unspecified domain.+ + * 'UFA.ERUnitFnApprox' (extends 'ERFnApprox'): generalises function graph enclosures+ on the domain @[-1,1]^n@.+ (/Not exported here, used only internally./)+ + * 'ERFnDomApprox' (extends 'ERFnApprox'):+ generalises function enclosures over a specified and queriable domain box+ (instance of class 'DomainBox').+ + At all levels, all field operations are supported as well as+ some elementary operations, namely exp, sin and cos.+ Log and sqrt are planned to be added soon. + + Implementations of 'UFB.ERUnitFnBase':+ + * 'ERChebPoly' + + By using the Chebyshev basis on domain @[-1,1]^n@, + we gain simple and optimally rounding degree reduction + as well as relatively simple handling of rounding + in other operations.++ Implementations of 'UFA.ERUnitFnApprox':+ + * 'ERFnInterval'+ + Implementations of 'ERFnDomApprox':+ + * 'ERFnDomTranslApprox': + builds a basic implementation + using an instance of 'UFA.ERUnitFnApprox'.+ + * 'ERFnTuple': + extends another implementation of 'ERFnDomApprox' + to work with lists of functions simultaneously.+ + * 'ERFnDomEdgesApprox': + separately enclose a function on its domain box + as well as on all the domain's hyper-edges + (including the corners) using+ another implementation of 'ERFnDomApprox'.+ + * 'ERFnPiecewise': + allows the domain box to be bisected + to an arbitrary finite depth + and uses another implementation of 'ERFnDomApprox' + to approximate the function on each segment. +-}+module Data.Number.ER.RnToRm +(+ module Data.Number.ER.RnToRm.DefaultRepr,+ module Data.Number.ER.RnToRm.Approx,+ module Data.Number.ER.Real.DomainBox+)+where++import Data.Number.ER.RnToRm.DefaultRepr+import Data.Number.ER.RnToRm.Approx+import Data.Number.ER.Real.DomainBox++import qualified Data.Number.ER.RnToRm.UnitDom.Approx as UFA+import qualified Data.Number.ER.RnToRm.UnitDom.Base as UFB++import qualified Data.Number.ER.Real.Approx as RA+import qualified Data.Number.ER.Real.Approx.Elementary as RAEL++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom+import Data.Number.ER.RnToRm.UnitDom.Approx.Interval+import Data.Number.ER.RnToRm.Approx.DomTransl+import Data.Number.ER.RnToRm.Approx.DomEdges+import Data.Number.ER.RnToRm.Approx.Tuple+import Data.Number.ER.RnToRm.Approx.PieceWise+
+ src/Data/Number/ER/RnToRm/Approx.hs view
@@ -0,0 +1,303 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-|+ Module : Data.Number.ER.RnToRm.Approx+ Description : classes abstracting function approximations+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++ Approximation of a real functions with rectangular domains.+ + To be imported qualified, usually with the synonym FA. +-}+module Data.Number.ER.RnToRm.Approx+(+ ERFnApprox(..),+ ERFnDomApprox(..),+ bisectUnbisectDepth+)+where++import Prelude hiding (const)++import qualified Data.Number.ER.Real.Approx as RA+import qualified Data.Number.ER.Real.DomainBox as DBox+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)+import Data.Number.ER.BasicTypes++import qualified Data.Map as Map++{-|+ A class of types that approximate first-order real functions+ @R^n -> R^m@ using some type of graph enclosures. The domains+ of the functions can be neither specified nor investigated + by operations in this class.++ This class extends 'RA.ERApprox' so that we could perform point-wise+ operations on the functions.++ This class is associated with:+ + * two real number types (instances of 'RA.ERIntApprox') for working with parts of the function's domain and range;+ + * a type of boxes indexed by variables (instance of 'DomainBox') for working with+ parts of the function's domain.+-}+class + (RA.ERApprox fa, RA.ERIntApprox domra, RA.ERIntApprox ranra, + DomainBox box varid domra) => + ERFnApprox box varid domra ranra fa+ | fa -> box varid domra ranra+ where+ {-| Check internal consistency and report problem if any. -}+ check :: + String {-^ indentification of caller location for easier debugging -} -> + fa -> fa+ domra2ranra :: + fa {-^ this parameter is not used except for type checking -} -> + domra -> ranra+ ranra2domra :: + fa {-^ this parameter is not used except for type checking -} -> + ranra -> domra + {-| + Get the internal degree of quality (usually polynomial degree) + of the approximation. + -}+ getDegree :: fa -> Int+ {-| + Set an upper bound on the degree of this function approximation.+ + This reduces the degree immediately if necessary and also+ affects all operations performed with this value later.+ -}+ setMaxDegree :: Int -> fa -> fa+ {-| + Get the current uppend bound on the degree associated + with this function approximation. + -}+ getMaxDegree :: fa -> Int+ {-| + Give a close upper bound of the precision of the range + at the best approximated point in the domain.+ -}+ getBestPrecision :: fa -> Precision+ {-| + Find some upper and lower bounds of the function over @[-1,1]^n@.+ -} + getRangeApprox :: fa -> ranra+ {-| + Combine several functions with the same domain into one /tuple function/. + -}+ tuple :: [fa] -> fa+ {-|+ Reveal how many functions are bundled together.+ -}+ getTupleSize :: fa -> Int+ {-| + Modify a tuple of functions in a way + that does not treat the tuple elements uniformly.+ -}+ applyTupleFn ::+-- (ERFnApprox box varid domra ranra fa2) => +-- ([fa2] -> [fa2]) -> (fa -> fa)+ ([fa] -> [fa]) -> (fa -> fa)+ {-| + Find close upper and lower bounds of the volume of the entire enclosure.+ A negative volume means that the enclosure is certainly inconsistent.+ -} + volume :: fa -> ranra+ {-|+ Multiply a function approximation by a real number approximation.+ -}+ scale :: ranra -> fa -> fa+ {-|+ Intersect one enclosure by another but only on a box within its domain.+ -}+ partialIntersect ::+ EffortIndex -> + box {-^ the subdomain; defined by clipping the range of some variables -} ->+ fa {-^ function to improve by intersecting its subdomain -} -> + fa {-^ the enclosure to be used on the subdomain (but defined on the whole domain) -} ->+ fa+ {-|+ Intersect two enclosures and measure the global improvement as one number.+ + (Use 'RA.intersectMeasureImprovement' defined in module "Data.Number.ER.Real.Approx" + to measure the improvement using a function enclosure.) + -} + intersectMeasureImprovement ::+ EffortIndex -> + fa -> + fa -> + (fa, ranra)+ {-^ enclosure intersection and measurement of improvement analogous to the one + returned by pointwise 'intersectMeasureImprovement' -}+ {-|+ Evaluate the function at the given point.+ -}+ eval :: box -> fa -> [ranra]+ {-|+ Fix some variables in the function to the given exact values.+ -}+ partialEval :: box -> fa -> fa+ {-| + A simple and limited composition of functions.+ + It is primarily intended to be used for precomposition with affine functions.+ -} + composeThin ::+ fa {-^ enclosure of @f@ -} ->+ Map.Map varid fa+ {-^ specifies the variables to substitute and for each such variable @v@, + gives an /exact/ enclosure of a function @f_v@ to substitute for @v@ -} ->+ fa + {-^ enclosure of @f[v |-> f_v]@ + + BEWARE: Enclosure is probably incorrect where values of @f_v@ are outside the domain of @v@ in @f@.+ -}+++{-|+ This class extends 'ERFnApprox' by:+ + * making the domain of the function enclosure available for inspection;+ + * allowing the construction of basic function enclosures+ where the domain has to be specified.+-}+class + (ERFnApprox box varid domra ranra fa,+ DomainIntBox box varid domra) => + ERFnDomApprox box varid domra ranra fa+ | fa -> box varid domra ranra+ where+ {-| + A function enclosure with no information about the function's values.+ -}+ bottomApprox :: + box {-^ the domain of the function -} -> + Int {-^ how many functions are bundled in this tuple -} -> + fa+ {-|+ Construct a constant enclosure for a tuple of functions.+ -}+ const :: box -> [ranra] -> fa+ {-|+ Construct the exact enclosure for a projection function+ (ie a variable).+ -}+ proj :: box -> varid -> fa+ {-|+ Return the domain of the function enclosure.+ -}+ dom :: fa -> box+ {-| + Split the domain into two halves, yoelding two function enclosures.+ -}+ bisect :: + varid {-^ variable (axis) to split on -} -> + Maybe domra {-^ where exactly to split (this has to be exact) -} -> + fa -> + (fa, fa)+ {-| + Merge function enclosures with neighbouring domains.+ -}+ unBisect :: + varid {-^ variable (axis) to glue on -} -> + (fa, fa) -> + fa+ {-| + Safely integrate a @R^n -> R^m@ function enclosure+ with some initial condition (origin and function at origin).+ -} + integrate :: + EffortIndex {-^ how hard to try -} ->+ fa {-^ function to integrate -} ->+ varid {-^ @x@ = variable to integrate by -} ->+ box {-^ integration range -} ->+ domra {-^ origin in terms of @x@; this has to be thin! -} ->+ fa {-^ values at origin -} ->+ fa+ {-| + Safely integrate a @R -> R^m@ function enclosure.+ -} + integrateUnary :: + EffortIndex {-^ how hard to try -} ->+ fa {-^ unary function to integrate -} ->+ domra {-^ integration range -} ->+ domra {-^ origin -} ->+ [ranra] {-^ values at origin -} ->+ fa+ -- default implementation reduces this to integrateMeasureImprovement:+ integrateUnary ix fD support origin vals =+ integrate ix fD defaultVar (DBox.unary support) origin (const (DBox.noinfo) vals)+ {-| + Safely integrate a @R^n -> R^m@ function enclosure+ intersecting it with a prior enclosure for the result.+ + The prior enclosure could contains one of more initial value.+ -} + integrateMeasureImprovement :: + EffortIndex {-^ how hard to try -} ->+ fa {-^ function to integrate -} ->+ varid {-^ variable to integrate by -} ->+ box {-^ integration domain -} ->+ domra + {-^ a sub-domain with relevant new information - + either about initial value(s) or about derivative -} ->+ fa {-^ approximation to result, including initial value(s) -} -> + (fa, fa) + {-^ improved result and measurement of improvement analogous to the one + returned by pointwise 'intersectMeasureImprovement' -}+ {-| + Safely integrate a @R -> R^m@ function enclosure+ intersecting it with a prior enclosure for the result.+ + The prior enclosure could contains one of more initial value.+ -} + integrateMeasureImprovementUnary :: + EffortIndex {-^ how hard to try -} ->+ fa {-^ unary function to integrate -} ->+ domra {-^ integration domain -} ->+ domra + {-^ a sub-domain with relevant new information - + either about initial value(s) or about derivative -} ->+ fa {-^ approximation to result, including initial value(s) -} -> + (fa, fa) + {-^ improved result and measurement of improvement analogous to the one + returned by pointwise 'intersectMeasureImprovement' -}+ -- default implementation reduces this to integrateMeasureImprovement:+ integrateMeasureImprovementUnary ix fD support origin fP =+ integrateMeasureImprovement ix fD defaultVar (DBox.unary support) origin fP+ + +{-|+ Recursively perform a number of bisections and then+ glue the bits back together. + + This way we can ensure that+ a piece-wise enclosure has a partition that goes+ to at least the given depth. +-}+bisectUnbisectDepth ::+ (ERFnDomApprox box varid domra ranra fa) =>+ Int {-^ required depth of bisection -} ->+ fa -> + fa+bisectUnbisectDepth depth f =+ aux splitVars depth f+ where+ splitVars = concat $ repeat $ DBox.keys $ dom f+ aux (var : restVars) depthsToGo f + | depthsToGo <= 0 = f+ | otherwise =+ unBisect var (fLDone, fRDone)+ where+ fLDone = aux restVars depthsToGoM1 fL+ fRDone = aux restVars depthsToGoM1 fR+ (fL, fR) = bisect var Nothing f+ depthsToGoM1 = depthsToGo - 1
+ src/Data/Number/ER/RnToRm/Approx/DomEdges.hs view
@@ -0,0 +1,487 @@+{-# OPTIONS_GHC -fno-warn-missing-methods #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-|+ Module : Data.Number.ER.RnToRm.Approx.DomEdges+ Description : separate approximations per domain-box hyper-edge+ Copyright : (c) Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++-}+module Data.Number.ER.RnToRm.Approx.DomEdges +(+ ERFnDomEdgesApprox(..)+)+where++import qualified Data.Number.ER.RnToRm.Approx as FA+import qualified Data.Number.ER.Real.Approx as RA+import qualified Data.Number.ER.Real.Approx.Elementary as RAEL++import qualified Data.Number.ER.Real.DomainBox as DBox+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox)+import Data.Number.ER.BasicTypes+import Data.Number.ER.Misc+import Data.Number.ER.PlusMinus++import Data.Typeable+import Data.Generics.Basics+import Data.Binary++import qualified Data.Map as Map+import qualified Data.Set as Set+import Data.List++{-|+ Use some function approximation type and for each domain box+ keep a structure of function approximations of this type indexed+ by the hyper-edge structure. For each hyper-edge of the domain,+ the approximation has this edge as its domain.+ + E.g. for a 2D square domain there are:+ + * one approximation for the whole square+ + * four 1D approximations, one for each edge+ + * eight 0D approximations, one for each endpoint of each edge + -}+data ERFnDomEdgesApprox varid fa =+ ERFnDomEdgesApprox+ {+ erfnMainVolume :: fa,+ erfnEdges :: Map.Map (varid, PlusMinus) (ERFnDomEdgesApprox varid fa)+ }+ deriving (Typeable,Data)++instance (Ord a, Binary a, Binary b) => Binary (ERFnDomEdgesApprox a b) where+ put (ERFnDomEdgesApprox a b) = put a >> put b+ get = get >>= \a -> get >>= \b -> return (ERFnDomEdgesApprox a b)++edgesLift1 ::+ (fa -> fa) ->+ (ERFnDomEdgesApprox varid fa) -> (ERFnDomEdgesApprox varid fa)+edgesLift1 op (ERFnDomEdgesApprox mainEncl edges) = + ERFnDomEdgesApprox (op mainEncl) (Map.map (edgesLift1 op) edges)+ +edgesLift2 ::+ (Ord varid) =>+ (fa -> fa -> fa) ->+ (ERFnDomEdgesApprox varid fa) -> (ERFnDomEdgesApprox varid fa) -> (ERFnDomEdgesApprox varid fa)+edgesLift2 op f1@(ERFnDomEdgesApprox mainEncl1 edges1) f2@(ERFnDomEdgesApprox mainEncl2 edges2) + | Map.keys edges1 == Map.keys edges2 =+ ERFnDomEdgesApprox (mainEncl1 `op` mainEncl2) $+ Map.intersectionWith (edgesLift2 op) edges1 edges2+ | otherwise =+ edgesLift2 op f1a f2a+ where+ (f1a, f2a) = unifyEdgeVariables f1 f2++unifyEdgeVariables ::+ (Ord varid) =>+ ERFnDomEdgesApprox varid fa ->+ ERFnDomEdgesApprox varid fa ->+ (ERFnDomEdgesApprox varid fa, ERFnDomEdgesApprox varid fa)+unifyEdgeVariables + f1@(ERFnDomEdgesApprox fa1 edges1) + f2@(ERFnDomEdgesApprox fa2 edges2) =+ (ERFnDomEdgesApprox fa1 edges1amended, + ERFnDomEdgesApprox fa2 edges2amended)+ where+ vars1 = Set.map fst $ Map.keysSet edges1+ vars2 = Set.map fst $ Map.keysSet edges2+ vars = Set.union vars1 vars2+ newVars1 = vars2 `Set.difference` vars1 + newVars2 = vars1 `Set.difference` vars2 + (ERFnDomEdgesApprox _ edges1amended) = + foldl (\f v -> addVarToEdges v f) f1 $ Set.toList newVars1+ (ERFnDomEdgesApprox _ edges2amended) = + foldl (\f v -> addVarToEdges v f) f2 $ Set.toList newVars2++addVarToEdges ::+ (Ord varid) =>+ varid ->+ ERFnDomEdgesApprox varid fa ->+ ERFnDomEdgesApprox varid fa +addVarToEdges var f@(ERFnDomEdgesApprox fa edges) =+ (ERFnDomEdgesApprox fa edgesNew)+ where+ edgesNew =+ Map.insert (var, Plus) f $ + Map.insert (var, Minus) f $ + Map.map (addVarToEdges var) edges+++instance + (FA.ERFnDomApprox box varid domra ranra fa, Ord varid, VariableID varid) =>+ Show (ERFnDomEdgesApprox varid fa)+ where+ show f@(ERFnDomEdgesApprox fa edges) =+ showAux [] f+ where+ showAux varAssignments (ERFnDomEdgesApprox fa edges) =+ edgeDescription +++ show fa +++ (concat $ map showEdge $ Map.toList edges)+ where+ edgeDescription + | null varAssignments =+ "\n>>>>> main enclosure: "+ | otherwise =+ "\n>>>>> edge" ++ showVarAssignments varAssignments ++ ": "+ showVarAssignments varAssignments =+ concat $ map showVarAssignment $ reverse varAssignments+ showVarAssignment (varID, val) =+ " " ++ showVar varID ++ "=" ++ show val+ showEdge ((varId, pm), faEdge) =+ showAux ((varId, varDomEndpoint) : varAssignments) faEdge+ where + varDomEndpoint =+ case pm of+ Minus -> varDomLo+ Plus -> varDomHi + (varDomLo, varDomHi) = RA.bounds varDom+ varDom = DBox.findWithDefault RA.bottomApprox varId domB+ domB = FA.dom fa++instance+ (FA.ERFnApprox box varid domra ranra fa) =>+ Eq (ERFnDomEdgesApprox varid fa)+ where+ (ERFnDomEdgesApprox fa1 edges1) == (ERFnDomEdgesApprox fa2 edges2) =+ fa1 == fa2++instance+ (FA.ERFnApprox box varid domra ranra fa, Ord fa) =>+ Ord (ERFnDomEdgesApprox varid fa)+ where+ compare (ERFnDomEdgesApprox fa1 edges1) (ERFnDomEdgesApprox fa2 edges2) =+ compare fa1 fa2++instance+ (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>+ Num (ERFnDomEdgesApprox varid fa)+ where+ fromInteger n = ERFnDomEdgesApprox (fromInteger n) Map.empty+ negate = edgesLift1 negate+ (+) = edgesLift2 (+)+ (*) = edgesLift2 (*)++instance + (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>+ Fractional (ERFnDomEdgesApprox varid fa)+ where+ fromRational r = ERFnDomEdgesApprox (fromRational r) Map.empty+ recip = edgesLift1 recip +++instance + (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>+ RA.ERApprox (ERFnDomEdgesApprox varid fa)+ where+ getGranularity (ERFnDomEdgesApprox mainEncl edges) =+ RA.getGranularity mainEncl+ setGranularity gran = edgesLift1 (RA.setGranularity gran) + setMinGranularity gran = edgesLift1 (RA.setMinGranularity gran)+ f1 /\ f2 = edgesLift2 (RA./\) f1 f2+ intersectMeasureImprovement ix + f1@(ERFnDomEdgesApprox mainEncl1 edges1) + f2@(ERFnDomEdgesApprox mainEncl2 edges2) + | Map.keys edges1 == Map.keys edges2 =+ (ERFnDomEdgesApprox mainEnclIsect edgesIsect,+ ERFnDomEdgesApprox mainEnclImpr edgesImpr)+ | otherwise =+ RA.intersectMeasureImprovement ix f1a f2a+ where+ (f1a, f2a) = unifyEdgeVariables f1 f2+ (mainEnclIsect, mainEnclImpr) =+ RA.intersectMeasureImprovement ix mainEncl1 mainEncl2+ edgesIsect = Map.map fst edgesIsectImpr+ edgesImpr = Map.map snd edgesIsectImpr+ edgesIsectImpr =+ Map.intersectionWith (RA.intersectMeasureImprovement ix) edges1 edges2 + leqReals fa1 fa2 =+ RA.leqReals (erfnMainVolume fa1) (erfnMainVolume fa2)++instance + (FA.ERFnDomApprox box varid domra ranra fa, RA.ERIntApprox fa, VariableID varid) =>+ RA.ERIntApprox (ERFnDomEdgesApprox varid fa)+ where+-- doubleBounds = :: ira -> (Double, Double) +-- floatBounds :: ira -> (Float, Float)+-- integerBounds :: ira -> (ExtendedInteger, ExtendedInteger)+ bisectDomain maybePt (ERFnDomEdgesApprox mainEncl edges) =+ (ERFnDomEdgesApprox mainEnclLo edgesLo,+ ERFnDomEdgesApprox mainEnclHi edgesHi)+ where+ (mainEnclLo, mainEnclHi) = RA.bisectDomain maybePtMainEncl mainEncl+ edgesLoHi = Map.intersectionWith RA.bisectDomain maybePtEdges edges+ edgesLo = Map.map fst edgesLoHi + edgesHi = Map.map snd edgesLoHi + (maybePtMainEncl, maybePtEdges) =+ case maybePt of+ Nothing -> + (Nothing, + Map.map (const Nothing) edges)+ Just (ERFnDomEdgesApprox mainEnclPt edgesPt) ->+ (Just mainEnclPt,+ Map.map Just edgesPt)+ bounds (ERFnDomEdgesApprox mainEncl edges) =+ (ERFnDomEdgesApprox mainEnclLo edgesLo,+ ERFnDomEdgesApprox mainEnclHi edgesHi)+ where+ (mainEnclLo, mainEnclHi) = RA.bounds mainEncl+ edgesLoHi = Map.map (RA.bounds) edges+ edgesLo = Map.map fst edgesLoHi + edgesHi = Map.map snd edgesLoHi+ f1 \/ f2 = edgesLift2 (RA.\/) f1 f2++instance + (FA.ERFnDomApprox box varid domra ranra fa, RAEL.ERApproxElementary fa, VariableID varid) =>+ RAEL.ERApproxElementary (ERFnDomEdgesApprox varid fa)+ where+ abs ix = edgesLift1 $ RAEL.abs ix+ exp ix = edgesLift1 $ RAEL.exp ix+ log ix = edgesLift1 $ RAEL.log ix+ sin ix = edgesLift1 $ RAEL.sin ix+ cos ix = edgesLift1 $ RAEL.cos ix+ +instance + (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>+ FA.ERFnApprox box varid domra ranra (ERFnDomEdgesApprox varid fa)+ where+ check prgLocation (ERFnDomEdgesApprox mainEncl edges) =+ ERFnDomEdgesApprox + (FA.check prgLocation mainEncl) + (Map.mapWithKey checkEdge edges)+ where+ checkEdge (var, pm) edgeFA =+ FA.check (prgLocation ++ showVar var ++ show pm ++ ": ") edgeFA+ domra2ranra fa d =+ FA.domra2ranra (erfnMainVolume fa) d+ ranra2domra fa r =+ FA.ranra2domra (erfnMainVolume fa) r+ setMaxDegree maxDegree = edgesLift1 (FA.setMaxDegree maxDegree)+ getTupleSize (ERFnDomEdgesApprox mainEncl _) =+ FA.getTupleSize mainEncl+ tuple [] = error "ERFnDomEdgesApprox: FA.tuple: empty list"+ tuple fs =+ foldl1 consFs fs + where+ consFs = edgesLift2 $ \a b -> FA.tuple [a,b]+ applyTupleFn tupleFn fn = (edgesLift1 $ FA.applyTupleFn tupleFnNoEdges) fn+ where+ tupleFnNoEdges fas =+ map erfnMainVolume $+ tupleFn $+ map (\fa -> ERFnDomEdgesApprox fa (makeEdges fa (erfnEdges fn))) + fas+ makeEdges fa oldEdges =+ Map.mapWithKey (makeVarPMEdge fa) oldEdges+ makeVarPMEdge fa (var, pm) oldEdge =+ ERFnDomEdgesApprox faNoVar $ makeEdges faNoVar (erfnEdges oldEdge)+ where+ faNoVar =+ FA.partialEval (DBox.singleton var domEndPt) fa+ domEndPt =+ case pm of Minus -> domL; Plus -> domR+ (domL, domR) = RA.bounds dom+ [dom] = DBox.elems $ FA.dom fa+ volume (ERFnDomEdgesApprox mainEncl edges) = FA.volume mainEncl+ scale ratio = edgesLift1 (FA.scale ratio)+ partialIntersect ix substitutions + f1@(ERFnDomEdgesApprox mainEncl1 edges1) + f2@(ERFnDomEdgesApprox mainEncl2 edges2) + | Map.keys edges1 == Map.keys edges2 =+ ERFnDomEdgesApprox (FA.partialIntersect ix substitutions mainEncl1 mainEncl2) $+ Map.intersectionWithKey partialIntersectEdge edges1 edges2+ | otherwise =+ FA.partialIntersect ix substitutions f1a f2a+ where+ (f1a, f2a) = unifyEdgeVariables f1 f2+ partialIntersectEdge (var, pm) edge1 edge2 + | withinSubstitutions =+ FA.partialIntersect ix substitutions edge1 edge2+ | otherwise = edge1+ where+ withinSubstitutions =+ (varDomEndpoint pm) `RA.refines` varVal+ where+ varVal =+ DBox.findWithDefault RA.bottomApprox var substitutions+ varDomEndpoint Minus = varDomLO+ varDomEndpoint Plus = varDomHI+ (varDomLO, varDomHI) = RA.bounds varDom+ varDom = DBox.lookup "DomEdges: partialIntersect: " var $ FA.dom mainEncl1 + eval ptBox (ERFnDomEdgesApprox mainEncl edges) + | null edgeVals =+ mainVal+ | otherwise =+ foldl1 (zipWith (RA./\)) edgeVals+ where+ mainVal = FA.eval ptBox mainEncl+ edgeVals = + concat $ map edgeEval $ Map.toList edges+ edgeEval ((x, sign), edgeFA) + | xPt `RA.refines` xDomLo && sign == Minus =+ [FA.eval ptBoxNoX edgeFA]+ | xPt `RA.refines` xDomHi && sign == Plus =+ [FA.eval ptBoxNoX edgeFA]+ | otherwise = []+ where+ (xDomLo, xDomHi) = RA.bounds xDom+ xDom = DBox.findWithDefault RA.bottomApprox x $ FA.dom mainEncl+ xPt = DBox.findWithDefault RA.bottomApprox x ptBox+ ptBoxNoX = DBox.delete x ptBox+ partialEval substitutions f@(ERFnDomEdgesApprox mainEncl edges) =+ (ERFnDomEdgesApprox mainEnclSubst edgesSubst)+ where+ mainEnclSubst = FA.partialEval substitutions mainEnclSelect+ edgesSubst = + Map.map (FA.partialEval substitutionsSelect) $+ Map.filterWithKey (\ (varID,_) _ -> varID `DBox.notMember` substitutionsSelect) edgesSelect+ (ERFnDomEdgesApprox mainEnclSelect edgesSelect, substitutionsSelect) = + foldl selectVar (f, substitutions) $ DBox.toList substitutions+ selectVar (fPrev@(ERFnDomEdgesApprox _ edgesPrev), substitutionsPrev) (varID, varVal)+ | varVal `RA.refines` varDomLo =+ (Map.findWithDefault fPrev (varID, Minus) edgesPrev, substitutionsNew) + | varVal `RA.refines` varDomHi =+ (Map.findWithDefault fPrev (varID, Plus) edgesPrev, substitutionsNew) + | otherwise = (fPrev, substitutionsPrev)+ where+ (varDomLo, varDomHi) = RA.bounds varDom+ varDom = DBox.findWithDefault RA.bottomApprox varID $ FA.dom mainEncl+ substitutionsNew = DBox.delete varID substitutionsPrev+ +instance + (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>+ FA.ERFnDomApprox box varid domra ranra (ERFnDomEdgesApprox varid fa)+ where+ dom (ERFnDomEdgesApprox mainEncl edges) = FA.dom mainEncl+ bottomApprox domB tupleSize =+ ERFnDomEdgesApprox (FA.bottomApprox domB tupleSize) $+ Map.fromList $ concat $+ map varEdges $ DBox.toList domB+ where + varEdges (varId, _) =+ [((varId, Minus), fEdge), ((varId, Plus), fEdge)]+ where+ fEdge = + FA.bottomApprox (DBox.delete varId domB) tupleSize+ const domB vals =+ ERFnDomEdgesApprox (FA.const domB vals) $+ Map.fromList $ concat $+ map varEdges $ DBox.toList domB+ where + varEdges (varId, _) =+ [((varId, Minus), fEdge), ((varId, Plus), fEdge)]+ where+ fEdge = + FA.const (DBox.delete varId domB) vals+ proj domB i =+ ERFnDomEdgesApprox mainEncl edges+-- Nothing ->+-- error $ +-- "DomEdges: projection index " ++ show i +-- ++ " out of range for domain " ++ show domB+ where+ mainEncl = FA.proj domB i+ edges =+ Map.fromList $ concat $ map makeVarEdges $ DBox.toList domB+ makeVarEdges (varID, varDom)+ | i == varID =+ [((varID, Minus), FA.const domNoVar [FA.domra2ranra mainEncl idomLo]),+ ((varID, Plus), FA.const domNoVar [FA.domra2ranra mainEncl idomHi])]+ | otherwise =+ [((varID, Minus), faNoVar),+ ((varID, Plus), faNoVar)]+ where+ domNoVar = DBox.delete varID domB+ (idomLo, idomHi) = RA.bounds idom+ idom = DBox.lookup "DomEdges: FA.proj: " i domB+ faNoVar = FA.proj domNoVar i+ bisect var maybePt f@(ERFnDomEdgesApprox mainEncl edges) + | varAbsent = (f,f)+ | otherwise =+ (ERFnDomEdgesApprox mainEnclLo edgesLo,+ ERFnDomEdgesApprox mainEnclHi edgesHi)+ where+ varAbsent =+ Map.notMember (var, Minus) edges+ (mainEnclLo, mainEnclHi) = FA.bisect var maybePt mainEncl+ pt = + case maybePt of + Nothing -> RA.defaultBisectPt varDom+ Just pt -> pt+ where+ varDom = + DBox.findWithDefault RA.bottomApprox var $ FA.dom mainEncl+ edgesLo =+ Map.insert (var, Minus) (edges Map.! (var, Minus)) $+ Map.insert (var, Plus) fAtPt $+ edgesLoNoVar+ edgesHi =+ Map.insert (var, Minus) fAtPt $+ Map.insert (var, Plus) (edges Map.! (var, Plus)) $+ edgesHiNoVar+ fAtPt = FA.partialEval (DBox.singleton var pt) f+ edgesLoNoVar = Map.map fst edgesLoHiNoVar+ edgesHiNoVar = Map.map snd edgesLoHiNoVar+ edgesLoHiNoVar = + Map.map (FA.bisect var maybePt) edgesNoVar+ edgesNoVar = + Map.delete (var, Plus) $ Map.delete (var, Minus) edges+ integrate ix fD x integdomBox origin fInit =+ ERFnDomEdgesApprox mainEncl edges+ where+ (ERFnDomEdgesApprox mainEnclD edgesD, + fInitWithX@(ERFnDomEdgesApprox _ edgesInitWithX)) = + unifyEdgeVariables fD fInit+ (ERFnDomEdgesApprox mainEnclInit edgesInit) = + Map.findWithDefault fInitWithX (x, Minus) edgesInitWithX + mainEncl = + FA.integrate ix mainEnclD x integdomBox origin mainEnclInit+ edges = + Map.insert (x, Minus) (FA.partialEval (DBox.singleton x xDomLo) fNoX) $ + Map.insert (x, Plus) (FA.partialEval (DBox.singleton x xDomHi) fNoX) $+ edgesNoX+ fNoX = ERFnDomEdgesApprox mainEncl edgesNoX+ edgesNoX =+ Map.intersectionWithKey integrEdge edgesD edgesInit+ (xDomLo, xDomHi) = RA.bounds xDom+ xDom = DBox.findWithDefault RA.bottomApprox x $ FA.dom fD+ integrEdge (varID, _) edgeD edgeInit =+ FA.integrate ix edgeD x (DBox.delete varID integdomBox) origin edgeInit+ + integrateMeasureImprovement ix fD x integdomBox xOrigin fP =+-- unsafePrint +-- ("DomEdges: integrateMeasureImprovement: faIntegrLo = " ++ show faIntegrLo) + (faIntegr, faImprovement)+ where+ faIntegr =+ faIntegrIsect+-- case RA.compareReals (FA.volume faIntegrIsect) (FA.volume faIntegrRaw) of+-- Just LT -> faIntegrIsect+-- _ -> faIntegrRaw -- this is wrong - forgets initial conditions!+ (faIntegrIsect, faImprovement) = + RA.intersectMeasureImprovement ix fP faIntegrRaw+ faIntegrRaw + | RA.isExact xOrigin = faIntegrLo+ | otherwise = faIntegrLo RA./\ faIntegrHi+ (xOriginLo, xOriginHi) = RA.bounds xOrigin+ faIntegrLo = + FA.integrate ix fD x integdomBox xOriginLo faPxLo + faPxLo = + FA.partialEval (DBox.singleton x xOriginLo) fP + faIntegrHi = + FA.integrate ix fD x integdomBox xOriginHi faPxHi + faPxHi = + FA.partialEval (DBox.singleton x xOriginHi) fP + +
+ src/Data/Number/ER/RnToRm/Approx/DomTransl.hs view
@@ -0,0 +1,495 @@+{-# LANGUAGE CPP #-}+{-# OPTIONS_GHC -fno-warn-missing-methods #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-|+ 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+ 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++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 + +
+ src/Data/Number/ER/RnToRm/Approx/PieceWise.hs view
@@ -0,0 +1,522 @@+{-# OPTIONS_GHC -fno-warn-missing-methods #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE DeriveDataTypeable #-}++{-|+ Module : Data.Number.ER.RnToRm.Approx.PieceWise+ Description : arbitrary precision piece-wise-something function enclosures+ Copyright : (c) Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++ Arbitrary precision piece-wise something + (eg linear, polynomial, rational) enclosures + of functions @R^n->R^m@.+ + The type of approximation within segments is specified+ by an instance of 'FA.ERFnDomApprox'.+ + The piece-wise construction defines another instance of 'FA.ERFnDomApprox'.+-}+module Data.Number.ER.RnToRm.Approx.PieceWise +(+ ERFnPiecewise(..)+)+where++import qualified Data.Number.ER.RnToRm.BisectionTree as BISTR+import qualified Data.Number.ER.RnToRm.BisectionTree.Integration as BTINTEG++import qualified Data.Number.ER.RnToRm.Approx as FA +import qualified Data.Number.ER.Real.Approx as RA+import qualified Data.Number.ER.Real.Approx.Elementary as RAEL++import qualified Data.Number.ER.Real.DomainBox as DBox+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)+import Data.Number.ER.BasicTypes+import Data.Number.ER.Misc++import Data.Typeable+import Data.Generics.Basics+import Data.Binary++import Data.Maybe++{-|+ Arbitrary precision piece-wise something + (eg linear, polynomial, rational) enclosures + of functions @R^n->R^m@.+ + The type of approximation within segments is specified+ by an instance of 'FA.ERFnDomApprox'.+ + The piece-wise construction defines another instance of 'FA.ERFnDomApprox'.+-}+data ERFnPiecewise box varid domra fa = + ERFnPiecewise (BISTR.BisectionTree box varid domra fa)+ deriving (Typeable, Data)+ +instance (Binary a, Binary b, Binary c, Binary d) => Binary (ERFnPiecewise a b c d) where+ put (ERFnPiecewise a) = put a+ get = get >>= \a -> return (ERFnPiecewise a)+ +pwLift1 ::+ (DomainBox box varid domra) =>+ (fa -> fa) ->+ (ERFnPiecewise box varid domra fa) -> + (ERFnPiecewise box varid domra fa)+pwLift1 op (ERFnPiecewise bistr) = + ERFnPiecewise (BISTR.mapWithDom (const op) bistr)+ +pwLift2 ::+ (RA.ERIntApprox domra, FA.ERFnDomApprox box varid domra ranra fa) =>+ (fa -> fa -> fa) ->+ EffortIndex ->+ (ERFnPiecewise box varid domra fa) -> + (ERFnPiecewise box varid domra fa) -> + (ERFnPiecewise box varid domra fa)+pwLift2 op ix f1@(ERFnPiecewise bistr1) f2@(ERFnPiecewise bistr2) =+ ERFnPiecewise $ + fromJust $ fst $ + BISTR.combineWith faSplit faSplit opBistr ix bistr1 bistr2+ where+ opBistr domB val1 val2 =+ (Just $ op val1 val2, [])+ +pwbistrZipWith ::+ (RA.ERIntApprox domra, FA.ERFnDomApprox box varid domra ranra fa) =>+ (fa -> fa -> res) ->+ EffortIndex ->+ (BISTR.BisectionTree box varid domra fa) ->+ (BISTR.BisectionTree box varid domra fa) ->+ (BISTR.BisectionTree box varid domra res)+pwbistrZipWith op ix bistr1 bistr2 =+ fromJust $ fst $ + BISTR.combineWith faSplit faSplit opBistr ix bistr1 bistr2 + where+ opBistr domB val1 val2 =+ (Just $ op val1 val2, [])++pwSplit ::+ (RA.ERIntApprox domra, DomainBox box varid domra) =>+ (fa -> (fa, fa)) ->+ (ERFnPiecewise box varid domra fa) -> (ERFnPiecewise box varid domra fa, ERFnPiecewise box varid domra fa)+pwSplit op f@(ERFnPiecewise bistr) = + (ERFnPiecewise bistr1, ERFnPiecewise bistr2)+ where+ bistr1 = BISTR.mapWithDom (const fst) bistr12 + bistr2 = BISTR.mapWithDom (const snd) bistr12 + bistr12 = BISTR.mapWithDom (const op) bistr++faSplit :: + (RA.ERIntApprox domra, FA.ERFnDomApprox box varid domra ranra fa) =>+ BISTR.ValueSplitter box varid domra fa+faSplit ix depth domB fa var pt = + FA.bisect var (Just pt) fa + +faCombine ::+ (RA.ERIntApprox domra, FA.ERFnDomApprox box varid domra ranra fa) =>+ BISTR.ValueCombiner box varid domra fa+faCombine ix depth (BISTR.Leaf _ _ v) = v+faCombine ix depth bistr =+ error "PieceWise: faCombine: not defined yet"+ +instance + (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>+ Show (ERFnPiecewise box varid domra fa)+ where+ show f@(ERFnPiecewise bistr) =+ "\nERFnPiecewise:" ++ show bistr++instance+ (FA.ERFnDomApprox box varid domra ranra fa) =>+ Eq (ERFnPiecewise box varid domra fa)+ where+ (ERFnPiecewise bistr1) == (ERFnPiecewise bistr2) =+ error $+ "ERFnPiecewise: Eq: not implemented yet"++instance+ (FA.ERFnDomApprox box varid domra ranra fa) =>+ Ord (ERFnPiecewise box varid domra fa)+ where+ compare (ERFnPiecewise bistr1) (ERFnPiecewise bistr2) =+ error $+ "ERFnPiecewise: Ord: not implemented yet"++instance+ (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>+ Num (ERFnPiecewise box varid domra fa)+ where+ fromInteger n = ERFnPiecewise $ BISTR.const DBox.noinfo (fromInteger n)+ negate = pwLift1 negate+ (+) = pwLift2 (+) 10+ (*) = pwLift2 (*) 10++instance+ (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>+ Fractional (ERFnPiecewise box varid domra fa)+ where+ fromRational r = ERFnPiecewise $ BISTR.const DBox.noinfo (fromRational r)+ recip = pwLift1 recip++instance + (FA.ERFnDomApprox box varid domra ranra fa, VariableID varid) =>+ RA.ERApprox (ERFnPiecewise box varid domra fa)+ where+ getGranularity (ERFnPiecewise bistr) =+ foldl max 10 $ map RA.getGranularity $ BISTR.collectValues bistr+ setGranularity gran = pwLift1 (RA.setGranularity gran) + setMinGranularity gran = pwLift1 (RA.setMinGranularity gran)+ f1 /\ f2 = pwLift2 (RA./\) 10 f1 f2+ intersectMeasureImprovement ix f1@(ERFnPiecewise bistr1) f2@(ERFnPiecewise bistr2) =+-- unsafePrint+-- (+-- "ERFnPiecewise: intersectMeasureImprovement:"+-- ++ "\n f1 = " ++ show f1+-- ++ "\n f2 = " ++ show f2+-- ++ "\n isect = " ++ show (ERFnPiecewise bistrIsect)+-- ++ "\n impr = " ++ show (ERFnPiecewise bistrImpr)+-- )+-- | length fas1 == length fas2 =+ (ERFnPiecewise bistrIsect, ERFnPiecewise bistrImpr)+-- | otherwise =+-- error $ show $ f1 RA./\ f2 + where+ bistrIsect = BISTR.mapWithDom (const fst) bistrIsectImpr + bistrImpr = BISTR.mapWithDom (const snd) bistrIsectImpr+ bistrIsectImpr = pwbistrZipWith (RA.intersectMeasureImprovement ix) ix bistr1 bistr2 + leqReals f1@(ERFnPiecewise bistr1) f2@(ERFnPiecewise bistr2) =+-- | length fas1 == length fas2 =+ leqTuple $ BISTR.collectValues $ pwbistrZipWith (RA.leqReals) 10 bistr1 bistr2+-- | otherwise =+-- error $ show $ f1 RA./\ f2+ where+ leqTuple [] = Just True+ leqTuple _ = + error $ "ERFnTuple: leqReals not implemented" ++instance + (FA.ERFnDomApprox box varid domra ranra fa, RA.ERIntApprox fa, VariableID varid) =>+ RA.ERIntApprox (ERFnPiecewise box varid domra fa)+ where+-- doubleBounds = :: ira -> (Double, Double) +-- floatBounds :: ira -> (Float, Float)+-- integerBounds :: ira -> (ExtendedInteger, ExtendedInteger)+ bisectDomain maybePt f@(ERFnPiecewise bistr) =+ case maybePt of+ Nothing ->+ pwSplit (RA.bisectDomain Nothing) f+ Just (ERFnPiecewise bistrPt) -> + (ERFnPiecewise bistr1, ERFnPiecewise bistr2)+ where+ bistr1 = BISTR.mapWithDom (const fst) bistr12 + bistr2 = BISTR.mapWithDom (const snd) bistr12 + bistr12 =+ pwbistrZipWith (\fa pt -> RA.bisectDomain (Just pt) fa) 10 + bistr bistrPt+ bounds = pwSplit RA.bounds+ f1 \/ f2 = pwLift2 (RA.\/) 10 f1 f2+ +instance + (FA.ERFnDomApprox box varid domra ranra fa, RAEL.ERApproxElementary fa, VariableID varid) =>+ RAEL.ERApproxElementary (ERFnPiecewise box varid domra fa)+ where+ abs ix = pwLift1 $ RAEL.abs ix+ exp ix = pwLift1 $ RAEL.exp ix+ log ix = pwLift1 $ RAEL.log ix+ sin ix = pwLift1 $ RAEL.sin ix+ cos ix = pwLift1 $ RAEL.cos ix+ +instance + (FA.ERFnDomApprox box varid domra ranra fa, + RA.ERIntApprox fa, + DomainBoxMappable box box varid domra domra, + Show box) =>+ FA.ERFnApprox box varid domra ranra (ERFnPiecewise box varid domra fa)+ where+ check prgLocation (ERFnPiecewise bistr) =+ ERFnPiecewise $ BISTR.mapWithDom checkSegm bistr+ where+ checkSegm dom f =+ FA.check (prgLocation ++ "segm " ++ show dom ++ ": ") f+ domra2ranra (ERFnPiecewise bistr) d =+ FA.domra2ranra fa d+ where+ (fa : _) = BISTR.collectValues bistr+ ranra2domra (ERFnPiecewise bistr) r =+ FA.ranra2domra fa r+ where+ (fa : _) = BISTR.collectValues bistr+ setMaxDegree maxDegree = pwLift1 (FA.setMaxDegree maxDegree)+ getTupleSize (ERFnPiecewise bistr) =+ FA.getTupleSize $ head $ BISTR.collectValues bistr+ tuple fs =+ foldl1 (pwLift2 (\a b -> FA.tuple [a,b]) 10) fs+ applyTupleFn tupleFn = pwLift1 $ FA.applyTupleFn tupleFnNoPW + where+ tupleFnNoPW fas =+ map (\ (ERFnPiecewise (BISTR.Leaf _ _ fa)) -> fa ) $+ tupleFn $+ map (\fa -> ERFnPiecewise $ BISTR.Leaf 0 (FA.dom fa) fa) + fas+ err = error "ERFnPiecewise: applyTupleFn"+ volume (ERFnPiecewise bistr) = + sum $ map FA.volume $ BISTR.collectValues bistr+ scale ratio = pwLift1 (FA.scale ratio)+ partialIntersect ix substitutions + f1@(ERFnPiecewise bistr1) + f2@(ERFnPiecewise bistr2) =+ ERFnPiecewise $ + head $+ BTINTEG.zipOnSubdomain + faSplit ix maxDepth substitutions+ updateInside updateTouch updateAway + [bistr1, bistr2]+ where+ maxDepth = effIx2int ix+ updateInside dom [val1, val2] =+ [FA.partialIntersect ix substitutions val1 val2]+ updateTouch = updateInside+ updateAway dom [val1, val2] =+ [val1]+ eval ptBox (ERFnPiecewise bistr) =+ foldl1 (zipWith (RA.\/)) $ + map (\fa -> FA.eval ptBox fa) $+ BISTR.collectValues $ BISTR.lookupSubtreeDom bistr ptBox + partialEval substitutions f@(ERFnPiecewise bistr) = + pwLift1 (FA.partialEval substitutions) (ERFnPiecewise bistrNoVars)+ where+ bistrNoVars =+ BISTR.removeVars substitutions bistr++instance + (FA.ERFnDomApprox box varid domra ranra fa, RA.ERIntApprox fa, Show box, + DomainBoxMappable box box varid domra domra) =>+ FA.ERFnDomApprox box varid domra ranra (ERFnPiecewise box varid domra fa)+ where+ dom (ERFnPiecewise bistr) = BISTR.bistrDom bistr+ bottomApprox domB tupleSize =+ ERFnPiecewise (BISTR.const domB $ FA.bottomApprox domB tupleSize)+ const domB vals =+ ERFnPiecewise $+ BISTR.const domB $ FA.const domB vals+ proj domB i =+ ERFnPiecewise $ BISTR.Leaf 0 domB $ FA.proj domB i + bisect var maybePt (ERFnPiecewise bistr) =+ (ERFnPiecewise bistrLo, ERFnPiecewise bistrHi)+ where+ (BISTR.Node _ _ _ _ bistrLo bistrHi) =+ BISTR.split faSplit 10 var pt DBox.noinfo bistr + pt =+ case maybePt of+ Nothing -> + RA.defaultBisectPt $ DBox.lookup "PieceWise: bisect: " var (BISTR.bistrDom bistr)+ Just pt -> pt+ unBisect var (ERFnPiecewise bistr1, ERFnPiecewise bistr2) =+ ERFnPiecewise $ + BISTR.Node (depth1 - 1) dom var domVarMid bistr1 bistr2+ where+ depth1 = BISTR.bistrDepth bistr1+ dom1 = BISTR.bistrDom bistr1+ dom2 = BISTR.bistrDom bistr2+ dom = DBox.unionWith (RA.\/) dom1 dom2+ domVarMid =+ snd $ RA.bounds $+ DBox.lookup "ERFnPiecewise: FA.unbisect: " var dom1+ integrate ix fD@(ERFnPiecewise bistrD) x integdomBox origin (ERFnPiecewise bistrInit) =+ ERFnPiecewise bistrIntegr+ where+ [bistrIntegr] =+ BTINTEG.zipFromOrigin+ faSplit faCombine+ ix x origin (Just integdomBox)+ zipOutsideRange+ shouldSplit+ integrateOriginHere+ integrateOriginLower+ integrateOriginHigher+ [bistrD, bistrInit]+ zipOutsideRange maybeFromL maybeFromR [bistrD, bistrInit] =+ unsafePrint+ (+ "ERFnPiecewise: integrateMeasureImprovement: zipOutsideRange: "+ ++ "\n domB = " ++ show domB+ ++ "\n bottomFn = " ++ show bottomFn+ )+ [bistrPadj]+ where+ (ERFnPiecewise bistrPadj) =+ case (maybeFromL, maybeFromR) of+ (Nothing, Nothing) -> bottomFn+ (Just faLO, Nothing) ->+ FA.partialIntersect ix + (DBox.singleton x domLO) + bottomFn + (ERFnPiecewise (BISTR.Leaf depth domB faLO)) + (Nothing, Just faHI) -> + FA.partialIntersect ix + (DBox.singleton x domHI) + bottomFn + (ERFnPiecewise (BISTR.Leaf depth domB faHI))+ bottomFn =+ ERFnPiecewise $ BISTR.Leaf depth domB $ FA.bottomApprox domB (FA.getTupleSize fD)+ (domLO, domHI) = + RA.bounds $ + DBox.lookup "ERFnPieceWise: integrate: zipOutsideRange: " x domB+ domB = BISTR.bistrDom bistrD+ depth = BISTR.bistrDepth bistrD+ shouldSplit ix depth _ _ _ =+ depth < (effIx2int ix)+ integrateOriginHere ix depth dom [faD, faInit] =+-- unsafePrint+-- (+-- "ERFnPiecewise: integrateMeasureImprovement: integrateOriginHere: "+-- ++ "\n dom = " ++ show dom+-- ++ "\n faLO = " ++ show faLO+-- ++ "\n faHI = " ++ show faHI+-- )+ (faLO, [faIntegr], faHI)+ where+ faIntegr = + FA.integrate ix faD x integdomBox origin faInit+ faLO =+ FA.partialEval (DBox.singleton x domLO) faIntegr+ faHI =+ FA.partialEval (DBox.singleton x domHI) faIntegr+ (domLO, domHI) = + RA.bounds $ + DBox.lookup "ERFnPieceWise: integrate: integrateOriginHere: " x dom+ integrateOriginLower ix depth dom faLO [faD, faInit] =+ ([faIntegr], faHI)+ where+ faIntegr = + FA.integrate ix faD x integdomBox domLO faLO+ faHI =+ FA.partialEval (DBox.singleton x domHI) faIntegr+ (domLO, domHI) = + RA.bounds $ + DBox.lookup "ERFnPieceWise: integrate: integrateOriginLower: " x dom+ integrateOriginHigher ix depth dom [faD, faInit] faHI =+ (faLO, [faIntegr])+ where+ faIntegr = + FA.integrate ix faD x integdomBox domHI faHI+ faLO =+ FA.partialEval (DBox.singleton x domLO) faIntegr+ (domLO, domHI) = + RA.bounds $ + DBox.lookup "ERFnPieceWise: integrate: integrateOriginHigher: " x dom++ integrateMeasureImprovement ix (ERFnPiecewise bistrD) x integdomBox origin (ERFnPiecewise bistrP) =+ (ERFnPiecewise bistrIsect, ERFnPiecewise bistrImpr)+ where+ [bistrIsect, bistrImpr] =+ BTINTEG.zipFromOrigin+ faSplit faCombine+ ix x origin (Just integdomBox)+ zipOutsideRange+ shouldSplit+ integrateOriginHere+ integrateOriginLower+ integrateOriginHigher+ [bistrD, bistrP]+ zipOutsideRange maybeFromL maybeFromR [bistrD, bistrP] =+-- unsafePrint+-- (+-- "ERFnPiecewise: zipOutsideRange"+-- )+ [bistrPadj, BISTR.mapWithDom (\d v -> FA.const d [1]) bistrP]+ where+ (ERFnPiecewise bistrPadj) =+ case (maybeFromL, maybeFromR) of+ (Nothing, Nothing) -> (ERFnPiecewise bistrP)+ (Just faLO, Nothing) ->+ FA.partialIntersect ix + (DBox.singleton x domLO) + (ERFnPiecewise bistrP) + (ERFnPiecewise (BISTR.Leaf depth domB faLO)) + (Nothing, Just faHI) -> + FA.partialIntersect ix + (DBox.singleton x domHI) + (ERFnPiecewise bistrP) + (ERFnPiecewise (BISTR.Leaf depth domB faHI)) + (domLO, domHI) = + RA.bounds $ + DBox.lookup "ERFnPieceWise: integrateMeasureImprovement: zipOutsideRange: " x domB+ domB = BISTR.bistrDom bistrP+ depth = BISTR.bistrDepth bistrP+ shouldSplit ix depth _ _ _ =+ depth < (effIx2int ix)+ integrateOriginHere ix depth dom [faD, faP] =+-- unsafePrint+-- (+-- "ERFnPiecewise: integrateMeasureImprovement: integrateOriginHere: "+-- ++ "\n dom = " ++ show dom+-- ++ "\n faLO = " ++ show faLO+-- ++ "\n faHI = " ++ show faHI+-- )+ (faLO, [faIsect, faImpr], faHI)+-- (FA.check "ERFnPieceWise: integrateOriginHere: faLO: " faLO, +-- [FA.check "ERFnPieceWise: integrateOriginHere: faIsect: " faIsect, +-- FA.check "ERFnPieceWise: integrateOriginHere: faImpr: " faImpr], +-- FA.check "ERFnPieceWise: integrateOriginHere: faHI: " faHI)+ where+ (faIsect, faImpr) = + FA.integrateMeasureImprovement ix faD x integdomBox origin faP+-- FA.integrateMeasureImprovement ix +-- (FA.check "ERFnPieceWise: integrateOriginHere: faD: " faD)+-- x integdomBox origin +-- (FA.check "ERFnPieceWise: integrateOriginHere: faP: " faP)+ faLO =+ FA.partialEval (DBox.singleton x domLO) faIsect+ faHI =+ FA.partialEval (DBox.singleton x domHI) faIsect+ (domLO, domHI) = + RA.bounds $ + DBox.lookup "ERFnPieceWise: integrateMeasureImprovement: integrateOriginHere: " x dom+ integrateOriginLower ix depth dom faLO [faD, faP] =+-- unsafePrint+-- (+-- "ERFnPiecewise: integrateMeasureImprovement: integrateOriginLower: "+-- ++ "\n dom = " ++ show dom+-- ++ "\n faLO = " ++ show faLO+-- ++ "\n faPadj = " ++ show faPadj+-- ++ "\n faHI = " ++ show faHI+-- )+ ([faIsect, faImpr], faHI)+ where+ (faIsect, faImpr) = + FA.integrateMeasureImprovement ix faD x integdomBox domLO faPadj+ faPadj =+ FA.partialIntersect ix (DBox.singleton x domLO) faP faLO+ faHI =+ FA.partialEval (DBox.singleton x domHI) faIsect+ (domLO, domHI) = + RA.bounds $ + DBox.lookup "ERFnPieceWise: integrateMeasureImprovement: integrateOriginLower: " x dom+ integrateOriginHigher ix depth dom [faD, faP] faHI =+-- unsafePrint+-- (+-- "ERFnPiecewise: integrateMeasureImprovement: integrateOriginHigher: "+-- ++ "\n dom = " ++ show dom+-- ++ "\n faLO = " ++ show faLO+-- ++ "\n faHI = " ++ show faHI+ -- )+ (faLO, [faIsect, faImpr])+ where+ (faIsect, faImpr) = + FA.integrateMeasureImprovement ix faD x integdomBox domHI faPadj+ faPadj =+ FA.partialIntersect ix (DBox.singleton x domHI) faP faHI+ faLO =+ FA.partialEval (DBox.singleton x domLO) faIsect+ (domLO, domHI) = + RA.bounds $ + DBox.lookup "ERFnPieceWise: integrateMeasureImprovement: integrateOriginHigher: " x dom+
+ src/Data/Number/ER/RnToRm/Approx/Tuple.hs view
@@ -0,0 +1,272 @@+{-# OPTIONS_GHC -fno-warn-missing-methods #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-|+ Module : Data.Number.ER.RnToRm.Approx.Tuples+ Description : a list of approximations over the same domain+ Copyright : (c) Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++ Lists of function approximations over the same domain. +-}+module Data.Number.ER.RnToRm.Approx.Tuple +(+ ERFnTuple(..)+)+where++import qualified Data.Number.ER.RnToRm.Approx as FA +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.BasicTypes++import Data.Typeable+import Data.Generics.Basics+import Data.Binary+++{-|+ A tuple of function approximations allowing one to get from + functions @R^n->R@ to a function @R^n -> R^m@.+-}+data ERFnTuple fa = + ERFnTuple { erfnTuple :: [fa] }+ deriving (Typeable, Data)+ +instance (Binary a) => Binary (ERFnTuple a) where+ put (ERFnTuple a) = put a+ get = get >>= \a -> return (ERFnTuple a)++tuplesLift1 ::+ (fa -> fa) ->+ (ERFnTuple fa) -> (ERFnTuple fa)+tuplesLift1 op (ERFnTuple fas) = + ERFnTuple (map op fas)+ +tuplesLift2 ::+ (Show fa) =>+ String ->+ (fa -> fa -> fa) ->+ (ERFnTuple fa) -> (ERFnTuple fa) -> (ERFnTuple fa)+tuplesLift2 callerLocation op f1@(ERFnTuple fas1) f2@(ERFnTuple fas2) + | length fas1 == length fas2 =+ ERFnTuple $ zipWith op fas1 fas2+ | otherwise =+ error $ + callerLocation ++ "incompatible lengths: " + ++ show (length fas1) ++ " != " ++ show (length fas2)+ ++ "\n first argument = \n" ++ show fas1+ ++ "\n second argument = \n" ++ show fas2++tuplesSplit ::+ (fa -> (fa, fa)) ->+ (ERFnTuple fa) -> (ERFnTuple fa, ERFnTuple fa)+tuplesSplit op f@(ERFnTuple fas) = + (ERFnTuple fas1, ERFnTuple fas2)+ where+ (fas1, fas2) = unzip $ map op fas++-- version with Map.Map:+--data ERFnTuple fa = +-- ERFnTuple (Map.Map varid fa)+-- deriving (Typeable, Data)+-- +--tuplesLift1 ::+-- (fa -> fa) ->+-- (ERFnTuple fa) -> (ERFnTuple fa)+--tuplesLift1 op (ERFnTuple fas) = +-- ERFnTuple (Map.map op fas)+-- +--tuplesLift2 ::+-- (fa -> fa -> fa) ->+-- (ERFnTuple fa) -> (ERFnTuple fa) -> (ERFnTuple fa)+--tuplesLift2 op f1@(ERFnTuple fas1) f2@(ERFnTuple fas2) +-- | Map.keys fas1 == Map.keys fas2 =+-- ERFnTuple $ Map.intersectionWith op fas1 fas2+-- | otherwise =+-- error $ +-- "ERFnTuple: incompatible keys: " +-- ++ show (Map.keys fas1) ++ "\n*****\n" ++ show (Map.keys fas2)+--+--tuplesSplit ::+-- (fa -> (fa, fa)) ->+-- (ERFnTuple fa) -> (ERFnTuple fa, ERFnTuple fa)+--tuplesSplit op f@(ERFnTuple fas) = +-- (ERFnTuple fas1, ERFnTuple fas2)+-- where+-- fas1 = Map.map fst fas12 +-- fas2 = Map.map snd fas12+-- fas12 = Map.map op fas+ + +instance + (FA.ERFnDomApprox box varid domra ranra fa) =>+ Show (ERFnTuple fa)+ where+ show f@(ERFnTuple fas) =+ concat $ map showFA $ zip [0,1..] fas+ where+ showFA (fnname, fa) =+ "\n>>> Function " ++ show fnname ++ ":" ++ show fa++instance+ (FA.ERFnApprox box varid domra ranra fa) =>+ Eq (ERFnTuple fa)+ where+ (ERFnTuple fas1) == (ERFnTuple fas2) =+ fas1 == fas2++instance+ (FA.ERFnApprox box varid domra ranra fa, Ord fa) =>+ Ord (ERFnTuple fa)+ where+ compare (ERFnTuple fas1) (ERFnTuple fas2) =+ compare fas1 fas2++instance+ (FA.ERFnDomApprox box varid domra ranra fa) =>+ Num (ERFnTuple fa)+ where+ fromInteger n = ERFnTuple [fromInteger n]+ negate = tuplesLift1 negate+ (+) = tuplesLift2 "ERFnTuple: +: " (+)+ (*) = tuplesLift2 "ERFnTuple: *: " (*)++instance + (FA.ERFnDomApprox box varid domra ranra fa) =>+ Fractional (ERFnTuple fa)+ where+ fromRational r = ERFnTuple [fromRational r]+ recip = tuplesLift1 recip +++instance + (FA.ERFnDomApprox box varid domra ranra fa) =>+ RA.ERApprox (ERFnTuple fa)+ where+ getGranularity (ERFnTuple fas) =+ foldl max 10 $ map RA.getGranularity fas+ setGranularity gran = tuplesLift1 (RA.setGranularity gran) + setMinGranularity gran = tuplesLift1 (RA.setMinGranularity gran)+ f1 /\ f2 = tuplesLift2 "ERFnTuple: /\\: " (RA./\) f1 f2+ intersectMeasureImprovement ix f1@(ERFnTuple fas1) f2@(ERFnTuple fas2)+ | length fas1 == length fas2 =+ (ERFnTuple fasIsect, ERFnTuple fasImpr)+ | otherwise =+ error $ show $ f1 RA./\ f2 + where+ (fasIsect, fasImpr) = unzip $ zipWith (RA.intersectMeasureImprovement ix) fas1 fas2 + leqReals f1@(ERFnTuple fas1) f2@(ERFnTuple fas2)+ | length fas1 == length fas2 =+ leqTuple $ zipWith RA.leqReals fas1 fas2+ | otherwise =+ error $ show $ f1 RA./\ f2+ where+ leqTuple [] = Just True+ leqTuple _ = + error $ "ERFnTuple: leqReals not implemented" ++instance + (FA.ERFnDomApprox box varid domra ranra fa, RA.ERIntApprox fa) =>+ RA.ERIntApprox (ERFnTuple fa)+ where+-- doubleBounds = :: ira -> (Double, Double) +-- floatBounds :: ira -> (Float, Float)+-- integerBounds :: ira -> (ExtendedInteger, ExtendedInteger)+ bisectDomain maybePt f@(ERFnTuple fas) =+ case maybePt of+ Nothing ->+ tuplesSplit (RA.bisectDomain Nothing) f+ Just (ERFnTuple fasPt) -> + (ERFnTuple fas1, ERFnTuple fas2)+ where+ (fas1, fas2) = + unzip $ + map (\(fa, pt) -> RA.bisectDomain (Just pt) fa) $ + zip fas fasPt+ bounds = tuplesSplit RA.bounds+ f1 \/ f2 = tuplesLift2 "ERFnTuple: \\/: " (RA.\/) f1 f2++instance + (FA.ERFnDomApprox box varid domra ranra fa, RAEL.ERApproxElementary fa) =>+ RAEL.ERApproxElementary (ERFnTuple fa)+ where+ abs ix = tuplesLift1 $ RAEL.abs ix+ exp ix = tuplesLift1 $ RAEL.exp ix+ log ix = tuplesLift1 $ RAEL.log ix+ sin ix = tuplesLift1 $ RAEL.sin ix+ cos ix = tuplesLift1 $ RAEL.cos ix+ +instance + (FA.ERFnDomApprox box varid domra ranra fa) =>+ FA.ERFnApprox box varid domra ranra (ERFnTuple fa)+ where+ check prgLocation (ERFnTuple fs) =+ ERFnTuple $ map checkComp $ zip [0..] fs+ where+ checkComp (n, f) =+ FA.check (prgLocation ++ "fn" ++ show n ++ ": ") f+ domra2ranra (ERFnTuple (fa:_)) d =+ FA.domra2ranra fa d+ ranra2domra (ERFnTuple (fa:_)) r =+ FA.ranra2domra fa r+ setMaxDegree maxDegree = tuplesLift1 (FA.setMaxDegree maxDegree)+ getTupleSize (ERFnTuple fas) = length fas+ tuple fs + | sameDomains doms = + ERFnTuple $ concat $ map erfnTuple fs+ | otherwise = + error $ + "ERFnTuple: FA.tuple: incompatible domains:\n " + ++ (unlines $ map show fs)+ where+ sameDomains [_] = True+ sameDomains (a : rest@(b : _)) =+ sameab && (sameDomains rest)+ where+ sameab =+ and $ map snd $ DBox.zipWithDefault RA.bottomApprox RA.equalApprox a b+ doms = map FA.dom fs+ applyTupleFn tupleFn (ERFnTuple fs) =+ FA.tuple $ tupleFn $ map (\fa -> ERFnTuple [fa]) fs+ volume (ERFnTuple fas) = sum $ map (FA.volume) fas+ scale ratio = tuplesLift1 (FA.scale ratio)+ partialIntersect ix substitutions =+ tuplesLift2 "ERFnTuple: partialIntersect: " $ FA.partialIntersect ix substitutions+ eval ptBox (ERFnTuple fas) =+ concat $ map (FA.eval ptBox) fas+ partialEval substitutions = tuplesLift1 $ FA.partialEval substitutions+ +instance + (FA.ERFnDomApprox box varid domra ranra fa) =>+ FA.ERFnDomApprox box varid domra ranra (ERFnTuple fa)+ where+ dom (ERFnTuple (fa:_)) = FA.dom fa+ bottomApprox domB tupleSize =+ ERFnTuple $ replicate tupleSize $ FA.bottomApprox domB 1+ const domB vals =+ ERFnTuple $ map (\v -> FA.const domB [v]) vals+ proj domB i =+ ERFnTuple [FA.proj domB i] ++ bisect var maybePt =+ tuplesSplit $ FA.bisect var maybePt+ integrate ix (ERFnTuple fasD) x integdomBox origin (ERFnTuple fasInit) =+ ERFnTuple $ map integ $ zip fasD fasInit+ where+ integ (faD, faInit) =+ FA.integrate ix faD x integdomBox origin faInit+ integrateMeasureImprovement ix (ERFnTuple fasD) x integdomBox origin (ERFnTuple fasP) =+ (ERFnTuple fasIsect, ERFnTuple fasImpr) + where+ (fasIsect, fasImpr) =+ unzip $ map integ $ zip fasD fasP+ integ (faD, faP) =+ FA.integrateMeasureImprovement ix faD x integdomBox origin faP
+ src/Data/Number/ER/RnToRm/BisectionTree.hs view
@@ -0,0 +1,665 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-|+ Module : Data.Number.ER.RnToRm.BisectionTree+ Description : hierarchical domain partitions + Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++ Defines a representation for recursive bisections of @R^n@+ by hyperplanes, each of which is perpendicular to a base axis.+ + Arbitrary data can be associated with the sections of a partition.+ + To be imported qualified, usually with prefix BISTR.+-}+module Data.Number.ER.RnToRm.BisectionTree +(+ BisectionTree(..),+ Depth,+ ValueSplitter,+ ValueCombiner,+ isLeaf,+ const,+ removeVars,+ sync2,+ syncMany,+ split,+ mapWithDom,+ mapLeaves,+ doBistr,+ doMap,+ doMapLeaves,+ combineWith,+ collectValues,+ collectDomValues,+ lookupSubtreeDom+)+where++import Prelude hiding (const, map, compare)+import qualified Prelude++import qualified Data.Number.ER.Real.Approx as RA+import qualified Data.Number.ER.Real.DomainBox as DBox+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)+import Data.Number.ER.BasicTypes++import Data.Number.ER.Misc++import Data.Typeable+import Data.Generics.Basics+import Data.Binary+--import BinaryDerive++import Data.Maybe+++{-|+ * The root of the tree often represents the whole @R^n@.+ + * Each node splits the parent's space into two using+ a specified variable (ie direction) and an optional splitting point.+ + * By default, a split is taken at the point defined by the method 'RA.bisect'.+-}+data BisectionTree box varid d v =+ Leaf+ {+ bistrDepth :: Depth,+ bistrDom :: box, -- ^ domain+ bistrVal :: v -- ^ value estimate+ }+ |+ Node+ {+ bistrDepth :: Depth, -- ^ depth of this node+ bistrDom :: box, -- ^ domain+ bistrDir :: varid, -- ^ direction to split in+ bistrPt :: d, -- ^ point that the split is at+ bistrLO :: BisectionTree box varid d v, -- ^ the half towards -Infty in split dir+ bistrHI :: BisectionTree box varid d v -- ^ the half towards +Infty in split dir+ }+ deriving (Typeable, Data)+ +type Depth = Int + +{- the following has been generated by BinaryDerive -}+instance (Binary a, Binary b, Binary c, Binary d) => Binary (BisectionTree a b c d) where+ put (Leaf a b c) = putWord8 0 >> put a >> put b >> put c+ put (Node a b c d e f) = putWord8 1 >> put a >> put b >> put c >> put d >> put e >> put f+ get = do+ tag_ <- getWord8+ case tag_ of+ 0 -> get >>= \a -> get >>= \b -> get >>= \c -> return (Leaf a b c)+ 1 -> get >>= \a -> get >>= \b -> get >>= \c -> get >>= \d -> get >>= \e -> get >>= \f -> return (Node a b c d e f)+ _ -> fail "no parse"+{- the above has been generated by BinaryDerive -}+ + +instance (VariableID varid, Show d, Show v, DomainBox box varid d) => + Show (BisectionTree box varid d v)+ where+ show = showBisectionTree show+ +showBisectionTree showValue =+ showB+ where+ showB (Leaf depth dom val) =+ "\n" +++ (concat (replicate (depth * 2) ".")) ++ "o "+ +++ (concatWith "," (Prelude.map showVD $ DBox.toList dom))+ +++ " |---> " ++ showValue val+ showB (Node depth dom dir pt lo hi) =+ "\n" +++ (concat (replicate (depth * 2) ".")) ++ "o "+ +++ (concatWith "," (Prelude.map showVD $ DBox.toList dom))+ +++ " //" ++ showVar dir ++ "\\\\"+ +++ (concat $ Prelude.map (showBisectionTree showValue) [lo,hi])+ showVD (v,d) =+ showVar v ++ "->" ++ show d++isLeaf ::+ BisectionTree box varid d v ->+ Bool +isLeaf (Leaf _ _ _) = True+isLeaf (Node _ _ _ _ _ _) = False++const ::+-- (DomainIntBox box varid d) =>+ box ->+ v ->+ BisectionTree box varid d v +const dom value =+ Leaf 0 dom value+ +{-| + value splitter function - parameters are: + depth, domain of value, value, variable to split by, + point to split at; returns the two split values+-} +type ValueSplitter box varid d v =+ (EffortIndex -> Depth -> box -> v -> varid -> d -> (v,v))++type ValueCombiner box varid d v = + (EffortIndex -> Depth -> (BisectionTree box varid d v) -> v) + + +split ::+ (RA.ERIntApprox d, DomainBox box varid d) =>+ ValueSplitter box varid d v ->+ EffortIndex ->+ varid {-^ variable @x@ to split by -} ->+ d {-^ point in domain of @x@ to split at -} ->+ box {-^ domain to lookup @x@ in if tree's domain does not have @x@ -} ->+ BisectionTree box varid d v ->+ BisectionTree box varid d v+split valSplitter ix splitDir splitPt fallbackDom bistr =+ resultBistr+ where+ resultBistr = spl bistr+ spl (Leaf depth dom val) =+ Node depth dom splitDir splitPt childLO childHI+ where+ childLO =+ Leaf depthInc domLO valLO+ childHI =+ Leaf depthInc domHI valHI + (valLO, valHI) = + valSplitter ix depth dom val splitDir splitPt+ depthInc = depth + 1+ domLO = + DBox.insert splitDir dirDomLO dom+ domHI = + DBox.insert splitDir dirDomHI dom+ (dirDomLO, dirDomHI) =+ RA.bisectDomain (Just splitPt) dirDom+ dirDom =+ DBox.findWithDefault + (DBox.lookup "BisectionTree: split: fallbackDom: " splitDir fallbackDom)+ splitDir dom+ spl bistr@(Node depth dom dir pt childLO childHI)+ | dir == splitDir =+ case RA.compareReals pt splitPt of+ Just LT -> -- split on lower half+ Node depth dom dir pt+ (Node depthInc domChildLO splitDir splitPt childLOsplitLO childLOsplitHI)+ childHI+ Just GT -> -- split on higher half+ Node depth dom dir pt+ childLO+ (Node depthInc domChildHI splitDir splitPt childHIsplitLO childHIsplitHI)+ _ -> bistr+ | otherwise = -- splitDir < dir =+ Node depth dom dir pt+ (Node + depthInc domChildLO splitDir splitPt childLOsplitLO childLOsplitHI)+ (Node + depthInc domChildHI splitDir splitPt childHIsplitLO childHIsplitHI)+ -- | dir < splitDir =+ -- Node depth dom dir childLOsplit childHIsplit+ where+ depthInc = depth + 1+ domChildLO = bistrDom childLO+ domChildHI = bistrDom childHI+ childLOsplit@(Node _ _ _ _ childLOsplitLO childLOsplitHI) = + spl childLO+ childHIsplit@(Node _ _ _ _ childHIsplitLO childHIsplitHI) = + spl childHI+ +{-|+ Apply a function to all values, thus creating a new tree.+-} +mapWithDom ::+ (DomainBox box varid d) =>+ (box -> v1 -> v2) ->+ BisectionTree box varid d v1 ->+ BisectionTree box varid d v2+mapWithDom f bistr@(Leaf _ dom val) =+ bistr { bistrVal = f dom val }+mapWithDom f bistr@(Node _ _ _ _ cLO cHI) =+ bistr + { + bistrLO = mapWithDom f cLO,+ bistrHI = mapWithDom f cHI+ }+ +{-|+ Apply a function to all values, thus creating a new tree.+-} +mapLeaves ::+ (BisectionTree box varid d v1 -> BisectionTree box varid d v2) ->+ BisectionTree box varid d v1 ->+ BisectionTree box varid d v2+mapLeaves f bistr@(Leaf _ dom val) =+ f bistr+mapLeaves f bistr@(Node _ _ _ _ cLO cHI) =+ bistr + { + bistrLO = mapLeaves f cLO,+ bistrHI = mapLeaves f cHI+ }+ +{-|+ Apply a function to all values, thus creating a list of new trees.+-} +mapMultiLeaves ::+ (BisectionTree box varid d v1 -> [BisectionTree box varid d v2]) ->+ BisectionTree box varid d v1 ->+ [BisectionTree box varid d v2]+mapMultiLeaves f bistr@(Leaf _ dom val) =+ f bistr+mapMultiLeaves f bistr@(Node _ _ _ _ cLO cHI) =+ Prelude.map (replaceChildren bistr) $ zip (mapMultiLeaves f cLO) (mapMultiLeaves f cHI)+ where+ replaceChildren bistr (newLO, newHI) =+ bistr + { + bistrLO = newLO,+ bistrHI = newHI+ }+ +{-|+ Perform a given action on all branches of a bisection tree, left to right.+ (optionally now going below the given depth)+-} +doBistr ::+ (box -> v -> IO ()) ->+ Maybe Int ->+ BisectionTree box varid d v ->+ IO ()+doBistr f Nothing bistr =+ m bistr+ where+ m (Node _ _ _ _ lo hi) =+ do+ m lo+ m hi+ m (Leaf _ dom val) =+ f dom val +doBistr f (Just maxDepth) bistr =+ m maxDepth bistr+ where+ m maxDepth (Node depth dom _ _ lo hi) + | maxDepth > 0 =+ do+ m (maxDepth - 1) lo+ m (maxDepth - 1) hi+ | otherwise =+ error $ "BisectionTree: doBistr: maxDepth (=" ++ show maxDepth ++ ") breached"+-- m err (Leaf depth dom val)+-- where+-- val = head $ collectValues lo+-- err = + m _ (Leaf _ dom val) =+ f dom val++{-|+ Perform a given action on all branches of a bisection tree, left to right.+ (optionally now going below the given depth)+-} +doMap ::+ (Depth -> box -> v -> IO v) ->+ Maybe Int ->+ BisectionTree box varid d v ->+ IO (BisectionTree box varid d v)+doMap f Nothing bistr =+ m bistr+ where+ m bistr@(Node _ _ _ _ lo hi) =+ do+ newLo <- m lo+ newHi <- m hi+ return $ bistr { bistrLO = newLo, bistrHI = newHi }+ m bistr@(Leaf depth dom val) =+ do+ newVal <- f depth dom val+ return $ bistr { bistrVal = newVal }+doMap f (Just maxDepth) bistr =+ m maxDepth bistr+ where+ m maxDepth bistr@(Node depth dom _ _ lo hi) + | maxDepth > 0 =+ do+ newLo <- m (maxDepth - 1) lo+ newHi <- m (maxDepth - 1) hi+ return $ bistr { bistrLO = newLo, bistrHI = newHi }+ | otherwise =+ error $ "BisectionTree: doBistr: maxDepth (=" ++ show maxDepth ++ ") breached"+-- m err (Leaf depth dom val)+-- where+-- val = head $ collectValues lo+-- err = + m _ bistr@(Leaf depth dom val) =+ do+ newVal <- f depth dom val+ return $ bistr { bistrVal = newVal }+ +{-|+ Perform a given action on all branches of a bisection tree, left to right+ with the option of further branching the tree.+ (optionally now going below the given depth)+-} +doMapLeaves ::+ (BisectionTree box varid d v -> IO (BisectionTree box varid d v)) ->+ Maybe Int ->+ BisectionTree box varid d v ->+ IO (BisectionTree box varid d v)+doMapLeaves f Nothing bistr =+ m bistr+ where+ m bistr@(Node _ _ _ _ lo hi) =+ do+ newLo <- m lo+ newHi <- m hi+ return $ bistr { bistrLO = newLo, bistrHI = newHi }+ m bistr@(Leaf depth dom val) =+ do+ f bistr+doMapLeaves f (Just maxDepth) bistr =+ m maxDepth bistr+ where+ m maxDepth bistr@(Node depth dom _ _ lo hi) + | maxDepth > 0 =+ do+ newLo <- m (maxDepth - 1) lo+ newHi <- m (maxDepth - 1) hi+ return $ bistr { bistrLO = newLo, bistrHI = newHi }+ | otherwise =+ error $ "BisectionTree: doBistr: maxDepth (=" ++ show maxDepth ++ ") breached"+-- m err (Leaf depth dom val)+-- where+-- val = head $ collectValues lo+-- err = + m _ bistr@(Leaf depth dom val) =+ do+ f bistr++removeVars ::+ (RA.ERIntApprox d, DomainIntBox box varid d, + DomainBoxMappable box box varid d d) =>+ box ->+ BisectionTree box varid d v -> + BisectionTree box varid d v +removeVars substitutions bistr =+ aux (bistrDepth bistr) bistr+ where+ aux depth (Leaf _ dom val) =+ Leaf depth domNoVars val+ where+ domNoVars =+ DBox.difference dom substitutions+ aux depth (Node _ dom v pt lo hi) + | v `DBox.member` substitutions =+ case (vVal `RA.refines` vDomLO, vVal `RA.refines` vDomHI) of+ (True, _) -> aux depth lo+ (_, True) -> aux depth hi+ | otherwise =+ Node depth domNoVars v pt loNoVars hiNoVars+ where+ vVal = DBox.lookup loc v substitutions+ vDomLO = DBox.lookup loc v $ bistrDom lo+ vDomHI = DBox.lookup loc v $ bistrDom hi+ loc = "RnToRm.BisectionTree: removeVars: "+ domNoVars =+ DBox.difference dom substitutions+ loNoVars = aux (depth + 1) lo + hiNoVars = aux (depth + 1) hi ++{-|+ Ensure both trees have equal structure at the top level:+ either they are all leaves or they all split at the same+ direction with the same splitting point.+ + Also, unify the domains at the top level.+-} +sync2 :: + (RA.ERIntApprox d, DomainIntBox box varid d) =>+ ValueSplitter box varid d v1 ->+ ValueSplitter box varid d v2 ->+ EffortIndex ->+ BisectionTree box varid d v1 -> + BisectionTree box varid d v2 -> + (BisectionTree box varid d v1, BisectionTree box varid d v2)+sync2 valSplitter1 valSplitter2 ix bistr1 bistr2 =+ case getPt bistr1 bistr2 of+ Nothing -> + unifyDom bistr1 bistr2+ Just (var, pt, dom) ->+ unifyDom+ (split valSplitter1 ix var pt dom bistr1)+ (split valSplitter2 ix var pt dom bistr2)+ where+ getPt bistr1 bistr2 + | isLeaf bistr1 && isLeaf bistr2 = Nothing+ | isLeaf bistr1 =+ Just (bistrDir bistr2, bistrPt bistr2, bistrDom bistr2) + | otherwise =+ Just (bistrDir bistr1, bistrPt bistr1, bistrDom bistr1) + unifyDom bistr1 bistr2 =+ (bistr1 { bistrDom = dom }, + bistr2 { bistrDom = dom })+ where+ dom =+ DBox.unify "RnToRm.BisectionTree: sync: " dom1 dom2+ dom1 = bistrDom bistr1 + dom2 = bistrDom bistr2 + +{-|+ Ensure all the trees have equal structure at the top level:+ either they are all leaves or they all split at the same+ direction with the same splitting point.+ + Also, unify the domains at the top level.+-} +syncMany :: + (RA.ERIntApprox d, DomainIntBox box varid d) =>+ ValueSplitter box varid d v ->+ EffortIndex ->+ [BisectionTree box varid d v] -> + [BisectionTree box varid d v]+syncMany valSplitter ix bistrs =+ case getPt bistrs of+ Nothing -> unifyDom bistrs+ Just (var, pt, dom) ->+ unifyDom $+ Prelude.map (split valSplitter ix var pt dom) bistrs+ where+ getPt [] = Nothing+ getPt (bistr : rest) + | isLeaf bistr = getPt rest+ | otherwise = Just (bistrDir bistr, bistrPt bistr, bistrDom bistr) + unifyDom bistrs =+ Prelude.map (setDom dom) bistrs+ where+ setDom dom bistr = bistr { bistrDom = dom }+ dom = + foldl (DBox.unify "RnToRm.BisectionTree: sync: ") DBox.noinfo $+ Prelude.map bistrDom bistrs + +{-|+ Combine two bisection trees using a given value combining function.+ Where necessary, leaves are split so that the resulting tree's structure+ is the union of the two argument tree structures. Such splitting of+ values in leaves is performed by the provided functions.+-}+combineWith ::+ (RA.ERIntApprox d, DomainIntBox box varid d) =>+ ValueSplitter box varid d v1+ {-^ value splitter function for tree 1 -} ->+ ValueSplitter box varid d v2+ {-^ value splitter function for tree 2 -} ->+ (box -> v1 -> v2 -> (Maybe res, aux)) + {-^ partial function to combine values with -} ->+ EffortIndex ->+ (BisectionTree box varid d v1) ->+ (BisectionTree box varid d v2) ->+ (Maybe (BisectionTree box varid d res), [aux])+combineWith valSplitter1 valSplitter2 f ix bistr1 bistr2 =+ combineAux bistr1sync bistr2sync+ where+ (bistr1sync, bistr2sync) = + sync2 valSplitter1 valSplitter2 ix bistr1 bistr2+ combineAux+ bistr1@(Leaf _ dom val1) + bistr2@(Leaf _ _ val2) =+ case f dom val1 val2 of+ (Nothing, aux) -> (Nothing, [aux])+ (Just val, aux) -> (Just $ bistr1 { bistrVal = val }, [aux])+ combineAux + bistr1@(Node _ dom _ _ lo1 hi1)+ bistr2@(Node _ _ _ _ lo2 hi2) =+ (+ Just $ bistr1 + {+ bistrLO = fromJust mbistrLO,+ bistrHI = fromJust mbistrHI+ }+ , + auxLO ++ auxHI+ )+ where+ (mbistrLO, auxLO) = combineAux lo1Sync lo2Sync+ (mbistrHI, auxHI) = combineAux hi1Sync hi2Sync+ (lo1Sync, lo2Sync) = + sync2 valSplitter1 valSplitter2 ix lo1 lo2+ (hi1Sync, hi2Sync) = + sync2 valSplitter1 valSplitter2 ix hi1 hi2+ +{-|+ return all values in leafs (except those within some CE subtree)+ as a list (from the leftmost to the rightmost)+-} +collectValues ::+ BisectionTree box varid b a -> [a]+collectValues (Leaf _ _ val) = [val]+collectValues (Node _ _ _ _ cLO cHI) =+ (collectValues cLO) ++ (collectValues cHI)+ +{-|+ return all values in leafs (except those within some CE subtree)+ as a list (from the leftmost to the rightmost)+-} +collectDomValues ::+ BisectionTree box varid d v -> [(box, v)]+collectDomValues (Leaf _ dom val) = [(dom,val)]+collectDomValues (Node _ _ _ _ cLO cHI) =+ (collectDomValues cLO) ++ (collectDomValues cHI)+ + +{-|+ linear ordering on bisection trees+-} +compare ::+ (Ord d, Ord varid) =>+ (v -> v -> Ordering) ->+ (BisectionTree box varid d v) ->+ (BisectionTree box varid d v) ->+ Ordering+compare compValues (Leaf _ _ _) (Node _ _ _ _ _ _) = LT+compare compValues (Node _ _ _ _ _ _) (Leaf _ _ _) = GT+compare compValues (Leaf _ _ val1) (Leaf _ _ val2) =+ compValues val1 val2+compare compValues + (Node _ _ dir1 pt1 lo1 hi1) + (Node _ _ dir2 pt2 lo2 hi2) =+ compareComposeMany $+ [Prelude.compare dir1 dir2,+ Prelude.compare pt1 pt2,+ compare compValues lo1 lo2,+ compare compValues hi1 hi2]++{-|+ lookup the smallest subtree whose domain covers the given rectangle+-}+lookupSubtreeDom ::+ (RA.ERIntApprox d, DomainBox box varid d) =>+ (BisectionTree box varid d v) ->+ box {-^ domain to look up within the tree -} ->+ (BisectionTree box varid d v)+lookupSubtreeDom origBistr dom = + lk origBistr+ where+ lk bistr@(Leaf _ _ _) = bistr+ lk bistr@(Node _ _ _ _ lo hi)+ | and $ Prelude.map snd $ DBox.zipWithDefault RA.bottomApprox (RA.refines) dom domHI = lk hi+ | and $ Prelude.map snd $ DBox.zipWithDefault RA.bottomApprox (RA.refines) dom domLO = lk lo+ | otherwise = bistr+ where+ domLO = bistrDom lo+ domHI = bistrDom hi++{-|+ Update a value on a given sub-domain,+ bisecting the tree if necessary to obtain+ a better fit for the domain, but not below+ a given depth limit.+ + With multiple domain dimensions, split the domain according to+ `DBox.bestSplit'.+-}+updateVal ::+ (RA.ERIntApprox d, DomainIntBox box varid d,+ DomainBoxMappable box box varid d d) =>+ ValueSplitter box varid d v ->+ EffortIndex ->+ Depth + {-^ depth limit -} ->+ box + {-^ domain to update on -} ->+ (box -> v -> v) + {-^ how to update values that intersect the above domain -} ->+ (BisectionTree box varid d v) ->+ (BisectionTree box varid d v)+updateVal valSplitter ix maxDepth updateDom updateFn bistr =+ upd bistr+ where+ upd bistr+ | noOverlap = bistr+ | edgeTouch && (isLeaf bistr) =+ updateLeaf bistr+ -- assuming we can update values on edges without+ -- influence on the interior+ | insideUpdateDom =+ mapLeaves updateLeaf bistr+ | depth >= maxDepth =+ mapLeaves updateLeaf bistr+ | otherwise = + -- divide and conquer:+ Node depth dom dir pt bistrLdone bistrRdone + where+ updateLeaf bistr =+ bistr { bistrVal = updateFn (bistrDom bistr) (bistrVal bistr) }+ noOverlap = + or $ Prelude.map RA.isEmpty $ DBox.elems domOverlap+ domOverlap = + DBox.intersectionWith (RA./\) dom updateDom+ insideUpdateDom = + and $ Prelude.map snd $ DBox.zipWith RA.refines dom updateDom+ edgeTouch =+ and $ Prelude.map snd $ DBox.zipWithDefaultSecond RA.bottomApprox endPointTouch dom updateDom+ endPointTouch i1 i2 =+ i1L == i2R || i1R == i2L+ where+ (==) = RA.eqSingletons+ (i1L, i1R) = RA.bounds i1+ (i2L, i2R) = RA.bounds i2 + depth = bistrDepth bistr+ dom = bistrDom bistr+ bistrLdone = upd bistrL+ bistrRdone = upd bistrR+ (Node _ _ _ _ bistrL bistrR) + | (isLeaf bistr) =+ split valSplitter ix dir pt DBox.noinfo bistr+ | otherwise = bistr + (dir, pt) =+ DBox.bestSplit dom+
+ src/Data/Number/ER/RnToRm/BisectionTree/Integration.hs view
@@ -0,0 +1,278 @@+{-|+ Module : Data.Number.ER.RnToRm.BisectionTree.Integration+ Description : abstract zipping of domain partitions used for integration+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + To be imported qualified, usually with prefix BTINTEG.+-}+module Data.Number.ER.RnToRm.BisectionTree.Integration +(+ zipFromOrigin, zipOnSubdomain+)+where++import qualified Data.Number.ER.RnToRm.BisectionTree as BISTR+import qualified Data.Number.ER.Real.Approx as RA++import qualified Data.Number.ER.Real.DomainBox as DBox+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)+import Data.Number.ER.BasicTypes+import Data.Number.ER.Misc++--import qualified Data.Sequence as Seq+--import qualified Data.Map as Map+import Data.Maybe++{-|+ Transform a bunch of bisection trees over the same domain + by "integrating" them in a very abstract sense. + The trees are unified in their splitting patterns in the process.+ By supplying certain parameters, this function can in fact+ perform numerical integration of piece-wise polynomial functions.+ + It can be also viewed as a "zipping+folding" operator over bisection trees that+ generates another bunch of bisection trees, synchronously traversing the original trees+ from a certain point on a selected axis outwards in both directions, + carrying some data along.+-}+zipFromOrigin ::+ (RA.ERIntApprox d, DomainIntBox box varid d, Show v1, Show v2, Show valPass) =>+ BISTR.ValueSplitter box varid d v1 ->+ BISTR.ValueCombiner box varid d v1 ->+ EffortIndex ->+ varid+ {-^ variable @x@ (ie axis or direction) to integrate in -} ->+ d + {-^ origin in terms of variable @x@ -} ->+ (Maybe (box))+ {-^ support, ie the domain on which to zip+ (automatically extended to include origin when projected to @x@) -} ->+ (Maybe valPass -> Maybe valPass -> [BISTR.BisectionTree box varid d v1] -> [BISTR.BisectionTree box varid d v2]) + {-^ what to do outside the support, + possibly being passed values from left/right+ when leaving the support -} ->+ (EffortIndex -> BISTR.Depth -> (box) -> [v1] -> [v2] -> Bool) + {-^ should a leaf be split? -} ->+ (EffortIndex -> BISTR.Depth -> (box) -> [v1] -> (valPass,[v2],valPass)) + {-^ integrator for a leaf containing the origin -} ->+ (EffortIndex -> BISTR.Depth -> (box) -> valPass -> [v1] -> ([v2], valPass))+ {-^ integrator over a leaf that sees the origin towards -infinity -} ->+ (EffortIndex -> BISTR.Depth -> (box) -> [v1] -> valPass -> (valPass, [v2])) + {-^ integrator over a leaf that sees the origin towards +infinity -} ->+ [BISTR.BisectionTree box varid d v1] + {-^ input functions -} ->+ [BISTR.BisectionTree box varid d v2]+ {-^ output functions+ + The number of output functions does not have to be + the same as the number of input functions. + -}+zipFromOrigin+ valSplitter valCombiner ix+ ivar origin maybeResultSupport outerValTransformer+ decideShouldSplit integrLeafOH integrLeafOL integrLeafOR+ bistrs =+ resultBistrs+ where+ (_, resultBistrs, _) = + integrateBistrOriginHere $ BISTR.syncMany valSplitter ix bistrs + maybeSupport = -- extend resultSupport to cover the origin+ fmap extendToOrigin maybeResultSupport+ where+ extendToOrigin domB =+ case DBox.member ivar domB of+ True -> DBox.insertWith (RA.\/) ivar origin domB+ False -> domB+ -- the following function is used when we know the origin is within the current sub-domain:+ integrateBistrOriginHere bistrs@((BISTR.Leaf depth dom _) : _)+ | decideShouldSplit ix depth dom vals integrVals = -- must descend+ integrateBistrOriginHere $ + map (BISTR.split valSplitter ix var pt dom) bistrs+ | otherwise =+ (Just lVal, map (\v -> BISTR.Leaf depth dom v) integrVals, Just rVal)+ where+ (var, pt) = DBox.bestSplit dom+ vals = map BISTR.bistrVal bistrs+ (lVal, integrVals, rVal) =+ integrLeafOH ix depth dom vals+ integrateBistrOriginHere bistrs@((BISTR.Node depth dom var pt lBounds rBounds):_)+ | origin `RA.refines` rDom =+-- unsafePrint +-- ("BTINTEG: integrateBistrOriginHere: rDom = " ++ show rDom ++ +-- " origin = " ++ show origin +++-- " lValHI = " ++ show lValHI +++-- " rValHI = " ++ show rValHI) + (lValHI, bistrsIntgHI, rValHI)+ | origin `RA.refines` lDom =+-- unsafePrint +-- ("BTINTEG: integrateBistrOriginHere: lDom = " ++ show lDom ++ +-- " origin = " ++ show origin +++-- " lValLO = " ++ show lValLO +++-- " rValLO = " ++ show rValLO)+ (lValLO, bistrsIntgLO, rValLO)+ | otherwise = -- origin overlaps both sides+ -- have to amalgamate these trees:+ integrateBistrOriginHere $+ map (\b -> BISTR.Leaf depth dom (valCombiner ix depth b)) bistrs+ where+ lDom = DBox.lookup "BTINTEG: zipFromOrigin: Here: L: " var (BISTR.bistrDom lBounds)+ rDom = DBox.lookup "BTINTEG: zipFromOrigin: Here: R: " var (BISTR.bistrDom rBounds)+ -- recursion when origin is entirely to the right of the centre:+ bistrsIntgHI = + zipWith + (\lo hi -> BISTR.Node depth dom var pt lo hi) + lBoundsIntgHI rBoundsIntgHI + (lValHIHI, rBoundsIntgHI, rValHI) =+ integrateBistrOriginHere $ + BISTR.syncMany valSplitter ix $ map BISTR.bistrHI bistrs+ (lValHI, lBoundsIntgHI) =+ integrateBistrOriginRight + (BISTR.syncMany valSplitter ix $ map BISTR.bistrLO bistrs) + lValHIHI+ -- recursion when origin is entirely to the left of the centre:+ bistrsIntgLO = + zipWith + (\lo hi -> BISTR.Node depth dom var pt lo hi) + lBoundsIntgLO rBoundsIntgLO + (lValLO, lBoundsIntgLO, rValLOLO) =+ integrateBistrOriginHere $ + BISTR.syncMany valSplitter ix $ map BISTR.bistrLO bistrs+ (rBoundsIntgLO, rValLO) =+ integrateBistrOriginLeft + rValLOLO + (BISTR.syncMany valSplitter ix $ map BISTR.bistrHI bistrs)+ -- the following function is used when we know + -- the origin is to the LEFT of the current sub-domain:+ integrateBistrOriginLeft Nothing bistrs = + -- previously detected as being outside the support+ (outerValTransformer Nothing Nothing bistrs, Nothing)+ integrateBistrOriginLeft (Just lVal) bistrs@(bistr:_)+ | (isJust maybeSupport) &&+ (and $ Prelude.map snd $ + DBox.zipWithDefaultSecond RA.bottomApprox RA.isInteriorDisjoint + (BISTR.bistrDom bistr) + (fromJust maybeSupport)) = + -- outside the integration domain + (outerValTransformer (Just lVal) Nothing bistrs, Nothing)+ integrateBistrOriginLeft (Just lVal) bistrs@((BISTR.Leaf depth dom _) : _)+ | decideShouldSplit ix depth dom vals integrVals = -- improve granularity by splitting+ integrateBistrOriginLeft (Just lVal) $ + map (BISTR.split valSplitter ix var pt dom) bistrs+ | otherwise = + (map (\v -> BISTR.Leaf depth dom v) integrVals, + Just rVal)+ where+ (var, pt) = DBox.bestSplit dom+ vals = map BISTR.bistrVal bistrs+ (integrVals, rVal) =+ integrLeafOL ix depth dom lVal vals+ integrateBistrOriginLeft mlVal bistrs@((BISTR.Node depth dom var pt _ _):_) =+ (bistrsIntg, mrVal2)+ where+ bistrsIntg = + zipWith (\lo hi -> BISTR.Node depth dom var pt lo hi) lBoundsINT rBoundsINT + (lBoundsINT, mrVal1) = + integrateBistrOriginLeft mlVal $ + BISTR.syncMany valSplitter ix $ map BISTR.bistrLO bistrs + (rBoundsINT, mrVal2) =+ integrateBistrOriginLeft mrVal1 $ + BISTR.syncMany valSplitter ix $ map BISTR.bistrHI bistrs +-- -- the following function is used when we know +-- -- the origin is to the RIGHT of the current sub-domain:+ integrateBistrOriginRight bistrs Nothing = + -- previously detected as being outside the support+ (Nothing, outerValTransformer Nothing Nothing bistrs)+ integrateBistrOriginRight bistrs@(bistr:_) (Just rVal)+ | (isJust maybeSupport) &&+ (and $ Prelude.map snd $ + DBox.zipWithDefaultSecond RA.bottomApprox RA.isInteriorDisjoint + (BISTR.bistrDom bistr) + (fromJust maybeSupport)) = + -- outside the integration domain + (Nothing, outerValTransformer Nothing (Just rVal) bistrs)+ integrateBistrOriginRight bistrs@((BISTR.Leaf depth dom _) : _) (Just rVal)+ | decideShouldSplit ix depth dom vals integrVals = -- improve granularity by splitting+ integrateBistrOriginRight + (map (BISTR.split valSplitter ix var pt dom) bistrs)+ (Just rVal)+ | otherwise = + (Just lVal,+ map (\v -> BISTR.Leaf depth dom v) integrVals)+ where+ (var, pt) = DBox.bestSplit dom+ vals = map BISTR.bistrVal bistrs+ (lVal, integrVals) =+ integrLeafOR ix depth dom vals rVal+ integrateBistrOriginRight bistrs@((BISTR.Node depth dom var pt _ _):_) mrVal =+ (mlVal2, bistrsIntg)+ where+ bistrsIntg = + zipWith (\lo hi -> BISTR.Node depth dom var pt lo hi) lBoundsINT rBoundsINT + (mlVal2, lBoundsINT) = + integrateBistrOriginRight + (BISTR.syncMany valSplitter ix $ map BISTR.bistrLO bistrs) mlVal1 + (mlVal1, rBoundsINT) =+ integrateBistrOriginRight + (BISTR.syncMany valSplitter ix $ map BISTR.bistrHI bistrs) mrVal ++{-|+ Zip a list of bisection trees in synchrony but do something+ else inside and not inside a given subdomain.+ + Further splitting at default points will be done up to the given depth+ in an attempt to separate the subdomain as well as possible.+ + If the subdomain is not properly isolated by the splitting at the+ maximum depth, splits are made at irregular points to ensure full isolation+ of the subdomain.+-}+zipOnSubdomain ::+ (RA.ERIntApprox d, DomainIntBox box varid d) =>+ BISTR.ValueSplitter box varid d v1 ->+ EffortIndex ->+ BISTR.Depth + {-^ depth limit -} ->+ box+ {-^ subdomain @sd@ -} ->+ (box -> [v1] -> [v2])+ {-^ what to do with values /inside/ @sd@ -} ->+ (box -> [v1] -> [v2])+ {-^ what to do with values /outside/ @sd@ but /touching/ it -} ->+ (box -> [v1] -> [v2])+ {-^ what to do with values /outside/ @sd@ -} ->+ [BISTR.BisectionTree box varid d v1] ->+ [BISTR.BisectionTree box varid d v2]+zipOnSubdomain valSplitter ix maxDepth sdom updateInside updateTouch updateAway bistrs =+ resultBistrs+ where+ resultBistrs = + zz $ BISTR.syncMany valSplitter ix bistrs+ zz bistrs@(BISTR.Leaf depth dom _ : _) + | intersect = + case depth < maxDepth of+ True ->+ zz $ map (BISTR.split valSplitter ix var pt dom) bistrs + False ->+ error "BTINTEG: zipOnSubdomain: maxDepth reached but irregular splitting not implemented yet"+ | away = lift updateAway+ | touch = lift updateTouch+ | inside = lift updateInside+ where+ (var, pt) = DBox.bestSplit dom+ lift updateFn =+ map (BISTR.Leaf depth dom) $ + updateFn dom $ + map BISTR.bistrVal bistrs+ (away, touch, intersect, inside) =+ DBox.classifyPosition dom sdom+ zz bistrs@(BISTR.Node depth dom var pt _ _ : _) =+ zipWith + (\bLO bHI -> BISTR.Node depth dom var pt bLO bHI) + (zz $ map BISTR.bistrLO bistrs) + (zz $ map BISTR.bistrHI bistrs) +
+ src/Data/Number/ER/RnToRm/BisectionTree/Path.hs view
@@ -0,0 +1,134 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-|+ Module : Data.Number.ER.RnToRm.BisectionTree.Path+ Description : addressing and modifying leaves+ Copyright : (c) Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++ Utilities for addressing and modifying leaves of + binary bisection trees.+-}+module Data.Number.ER.RnToRm.BisectionTree.Path where++import qualified Data.Number.ER.RnToRm.Approx as FA+import qualified Data.Number.ER.Real.Approx as RA+import Data.Number.ER.Real.DomainBox (VariableID(..))+import Data.Number.ER.BasicTypes++import Data.Typeable+import Data.Generics.Basics+import Data.Binary+--import BinaryDerive++{-|+ A path in a binary tree.+ It is used mainly in connection with 'BisectionTree.BisectionTree'.+-}+data BisecTreePath =+ BTP_H | BTP_R BisecTreePath | BTP_L BisecTreePath+ deriving (Eq, Typeable, Data)++{- the following has been generated by BinaryDerive -} +instance Binary BisecTreePath where+ put BTP_H = putWord8 0+ put (BTP_R a) = putWord8 1 >> put a+ put (BTP_L a) = putWord8 2 >> put a+ get = do+ tag_ <- getWord8+ case tag_ of+ 0 -> return BTP_H+ 1 -> get >>= \a -> return (BTP_R a)+ 2 -> get >>= \a -> return (BTP_L a)+ _ -> fail "no parse"+{- the above has been generated by BinaryDerive -} + +instance Show BisecTreePath + where+ show BTP_H = ""+ show (BTP_L rest) = "L" ++ show rest+ show (BTP_R rest) = "R" ++ show rest++instance Read BisecTreePath+ where+ readsPrec p ('L' : rest) =+ case readsPrec p rest of+ [(restParsed, s)] -> [(BTP_L restParsed, s)]+ _ -> []+ readsPrec p ('R' : rest) =+ case readsPrec p rest of+ [(restParsed, s)] -> [(BTP_R restParsed, s)]+ _ -> []+ readsPrec p s = [(BTP_H, s)]++{-|+ Assuming that bisection happens at default points as defined by+ 'RA.bisectDomain' and starts from the given root interval.+-}+path2dom ::+ (RA.ERIntApprox ira) =>+ ira {-^ root interval -} ->+ BisecTreePath -> + ira+path2dom rootdom path =+ p2d path rootdom+ where+ p2d BTP_H acc = acc+ p2d (BTP_L rest) acc = + p2d rest $ fst $ RA.bisectDomain Nothing $ acc+ p2d (BTP_R rest) acc = + p2d rest $ snd $ RA.bisectDomain Nothing $ acc+ +{-|+ A representation of a binary tree with a hole that+ can be efficiently filled.+-}+data FnZipper f+ = FnZ_H f+ | FnZ_L (FnZipper f) f+ | FnZ_R f (FnZipper f)++{-|+ Lookup a subdomain of a function according to a bisection path.+ Return the restrited function as well as a zipper that allows+ an efficient modification of the function on the looked up+ subdomain.+-} +lookupSubdomain ::+ (FA.ERFnDomApprox box varid domra ranra fa) =>+ fa ->+ BisecTreePath ->+ (fa, FnZipper fa)+lookupSubdomain fn BTP_H = (fn, FnZ_H fn)+lookupSubdomain fn (BTP_L restPath) =+ (resFn, FnZ_L subZipper hiFn)+ where+ (resFn, subZipper) = lookupSubdomain loFn restPath+ (loFn, hiFn) = FA.bisect defaultVar Nothing fn+lookupSubdomain fn (BTP_R restPath) =+ (resFn, FnZ_R loFn subZipper)+ where+ (resFn, subZipper) = lookupSubdomain hiFn restPath+ (loFn, hiFn) = FA.bisect defaultVar Nothing fn+ +{-|+ Modify a function in its subdomain as expressed by+ the zipper.+-}+updateFnZ ::+ (FA.ERFnDomApprox box varid domra ranra fa) =>+ (FnZipper fa) {-^ a function on a larger domain and a highlighted subdomain -} ->+ fa {-^ a function of the highlighted subdomain -} ->+ fa+updateFnZ (FnZ_H _) fn = fn+updateFnZ (FnZ_L loZipper hiFn) fn =+ FA.unBisect defaultVar (loFn, hiFn)+ where+ loFn = updateFnZ loZipper fn+updateFnZ (FnZ_R loFn hiZipper) fn =+ FA.unBisect defaultVar (loFn, hiFn)+ where+ hiFn = updateFnZ hiZipper fn
+ src/Data/Number/ER/RnToRm/DefaultRepr.hs view
@@ -0,0 +1,64 @@+{-|+ Module : Data.Number.ER.Real.DefaultRepr+ Description : concise names for default function representations+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : non-portable (requires fenv.h)++ This module supplies default instances for the real number and function classes+ described in "Data.Number.ER.RnToRm".+ + These classes form loosely coupled boundaries between abstraction layers.+ Nevertheless, we usually have particular implementations in mind, as shown here.+ + To preserve the intended loose coupling, please use these definitions+ only in functions that cannot infer from their input or output data which type of function enclosures+ they should use. Eg a function to add 1 to an enclosure should have the type:+ + > add1 :: (ERFnApprox box varid domra ranra fa) => fa -> fa+ > add1 f = f + 1+ + and /not/: @add1 :: FAPWP -> FAPWP@+ +-}+module Data.Number.ER.RnToRm.DefaultRepr+(+ module Data.Number.ER.RnToRm.DefaultRepr,+ module Data.Number.ER.Real.DomainBox.IntMap+)+where++import Data.Number.ER.Real.DefaultRepr+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.BasicTypes++import Data.Number.ER.Real.DomainBox.IntMap++import Data.Number.ER.RnToRm.UnitDom.Base+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom+import Data.Number.ER.RnToRm.UnitDom.Approx.Interval+import Data.Number.ER.RnToRm.Approx.DomTransl+import Data.Number.ER.RnToRm.Approx.DomEdges+import Data.Number.ER.RnToRm.Approx.Tuple+import Data.Number.ER.RnToRm.Approx.PieceWise++--import BinaryDerive++import qualified Data.Map as Map++type FAPU = ERFnInterval (ERChebPoly (Box Int) B) IRA +type FAPD = ERFnDomTranslApprox (Box (DomTransl IRA)) VarID FAPU IRA+type FAPT = ERFnTuple FAPD+type FAPE = ERFnDomEdgesApprox VarID FAPT+type FAPWP = ERFnPiecewise (Box IRA) VarID IRA FAPE++--type FA = FAPWL+type FA = FAPWP+
+ src/Data/Number/ER/RnToRm/TestingDefs.hs view
@@ -0,0 +1,72 @@+{-|+ Module : Data.Number.ER.Real.TestingDefs+ Description : definitions useful for testing in ghci+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : non-portable (requires fenv.h)++ A few definitions useful for testing the enclosures code, eg in ghci.+-}+module Data.Number.ER.RnToRm.TestingDefs where++import Data.Number.ER.RnToRm.DefaultRepr++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 qualified Data.Map as Map++fapuConst1 = (UFA.const [1]) :: FAPU++fapdConst1 = (FA.const DBox.noinfo [1]) :: FAPD+fapdConstU = (FA.const DBox.noinfo [(-1) RA.\/ 1]) :: FAPD+fapdConst01 = (FA.const DBox.noinfo [0 RA.\/ 1]) :: FAPD+fapd04X0 = (FA.proj (DBox.fromAscList [(0,0 RA.\/ 4)]) 0) :: FAPD+fapd13X0 = (FA.proj (DBox.fromAscList [(0,1 RA.\/ 3)]) 0) :: FAPD+fapd12X1 = (FA.proj (DBox.fromAscList [(1,1 RA.\/ 2)]) 1) :: FAPD+fapdUX0 = (FA.proj (DBox.fromAscList [(0,(-1) RA.\/ 1)]) 0) :: FAPD+fapdUX1 = (FA.proj (DBox.fromAscList [(1,(-1) RA.\/ 1)]) 1) :: FAPD++fapeConst1 = (FA.const DBox.noinfo [1]) :: FAPE+fapeConstU = (FA.const DBox.noinfo [(-1) RA.\/ 1]) :: FAPE+fapeConst01 = (FA.const DBox.noinfo [0 RA.\/ 1]) :: FAPE++fape13X0 = (FA.proj (DBox.fromAscList [(0,1 RA.\/ 3)]) 0) :: FAPE+fape12X1 = (FA.proj (DBox.fromAscList [(1,1 RA.\/ 2)]) 1) :: FAPE +fapeUX0 = (FA.proj (DBox.fromAscList [(0,(-1) RA.\/ 1)]) 0) :: FAPE+fapeUX1 = (FA.proj (DBox.fromAscList [(1,(-1) RA.\/ 1)]) 1) :: FAPE++fapeTestMult = (fapeUX0 + (FA.setMaxDegree 3 fapeConst01)) * (fapeConstU)+fapeMultiVar = (fapeUX0 + fapeUX1 * fapeUX0 + fapeUX1 * fapeUX1)+fapeTestPEval = FA.partialEval (DBox.fromList [(1,2 RA.\/ 3)]) fapeMultiVar++fapeUConst1 = (FA.const (DBox.unary $ (0)RA.\/1) [1]) :: FAPE+fapeUConst13 = (FA.const (DBox.unary $ (0)RA.\/1) [1 RA.\/ 3]) :: FAPE+fapeUConst13InitPt = FA.partialIntersect 1 (DBox.unary 0) fapeUConst13 fapeUConst1+ +fapwUUX0 = (FA.proj (DBox.fromAscList [(0,(1) RA.\/ 1)]) 0) :: FAPWP+fapwUUX1 = (FA.proj (DBox.fromAscList [(1,(-1) RA.\/ 1)]) 1) :: FAPWP++fapwUX0 = (FA.proj (DBox.fromAscList [(0,(0) RA.\/ 1)]) 0) :: FAPWP+fapwUX1 = (FA.proj (DBox.fromAscList [(1,(0) RA.\/ 1)]) 1) :: FAPWP++fapwUConst1 = (FA.const (DBox.noinfo) [1]) :: FAPWP+fapwUConst13 = (FA.const (DBox.unary $ (0)RA.\/1) [1 RA.\/ 3]) :: FAPWP+fapwUConst13InitPt = FA.partialIntersect 1 (DBox.unary 0) fapwUConst13 fapwUConst1 ++testIntegrE = + FA.integrateMeasureImprovement 1 (FA.setMaxDegree 0 fapeUConst13InitPt) 0 (DBox.noinfo) 0 fapeUConst13InitPt++testIntegrP = + FA.integrateMeasureImprovement 1 (FA.setMaxDegree 0 fapwUConst13InitPt) 0 (DBox.unary $ 0 RA.\/ 0.5) 0 fapwUConst13InitPt++x = FA.setMaxDegree 4 fapwUX0+fn1 = (1 + x) RA.\/ (1 + 3*x)+fn2 = FA.integrateUnary 0 fn1 0 (0 RA.\/ 1) [1]+fn3 = FA.integrateUnary 0 fn2 0 (0 RA.\/ 1) [1] -- this seems wrong!
+ src/Data/Number/ER/RnToRm/UnitDom/Approx.hs view
@@ -0,0 +1,92 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.Approx+ Description : class abstracting function enclosures on @[-1,1]^n@+ Copyright : (c) Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++ Approximation of continuous real functions + defined on the unit rectangle domain of a certain dimension.+ + To be imported qualified, usually with the synonym UFA. +-}+module Data.Number.ER.RnToRm.UnitDom.Approx+(+ ERUnitFnApprox(..)+)+where++import Data.Number.ER.RnToRm.Approx+import qualified Data.Number.ER.Real.DomainBox as DBox+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)+import Data.Number.ER.BasicTypes++import qualified Data.Map as Map++{-|+ This class extends 'ERFnApprox' by:+ + * assuming that the domain of the function enclosures is always @[-1,1]^n@ for some @n@;+ + * allowing the construction of basic function enclosures+ where the domain has to be known.+-}++class (ERFnApprox box varid domra ranra fa) => + ERUnitFnApprox box varid domra ranra fa+ | fa -> box varid domra ranra+ where+ {-| + A function enclosure with no information about the function's values.+ -}+ bottomApprox :: fa+ {-|+ Construct a constant enclosure for a tuple of functions.+ -}+ const :: [ranra] -> fa+ {-| + Construct the exact enclosure of an affine function on @[-1,1]^n@. + -} + affine :: + [ranra] {-^ values at 0 -} ->+ Map.Map varid ([ranra]) {-^ ascents of each base vector -} -> + fa+ {-| + Find close upper and lower bounds of the volume of the entire enclosure.+ A negative volume means that the enclosure is certainly inconsistent.+ + Explicitly specify the variables to identify the dimension of the domain.+ -} + volume :: [varid] -> fa -> ranra+ {-|+ Intersect two enclosures and measure the global improvement as one number.+ + (Use 'RA.intersectMeasureImprovement' defined in module "Data.Number.ER.Real.Approx" + to measure the improvement using a function enclosure.) + + Explicitly specify the variables to identify the dimension of the domain.+ -} + intersectMeasureImprovement ::+ EffortIndex -> + [varid] ->+ fa -> + fa -> + (fa, ranra)+ {-^ enclosure intersection and measurement of improvement analogous to the one + returned by the pointwise 'RA.intersectMeasureImprovement' -}+ {-| + Safely integrate a @[-1,1]^n -> R^m@ function enclosure+ with some initial condition (origin and function at origin).+ -} + integrate :: + EffortIndex {-^ how hard to try -} ->+ fa {-^ function to integrate -} ->+ varid {-^ @x@ = variable to integrate by -} ->+ domra {-^ origin in terms of @x@; this has to be exact! -} ->+ fa {-^ values at origin -} ->+ fa
+ src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs view
@@ -0,0 +1,589 @@+{-# OPTIONS_GHC -fno-warn-missing-methods #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE DeriveDataTypeable #-}++{-|+ Module : Data.Number.ER.RnToRm.UnitDom.Approx.Interval+ Description : arbitrary precision function enclosures on @[-1,1]^n@+ Copyright : (c) Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++ A construction of an enclosure of a real function on+ the domain [-1,1]^n for some n using elements of some+ base (eg rational functions or polynomials).+-}+module Data.Number.ER.RnToRm.UnitDom.Approx.Interval +(+ ERFnInterval(..),+ ERFnContext(..)+)+where++import qualified Data.Number.ER.Real.Base as B+import Data.Number.ER.Real.Approx.Interval+import Data.Number.ER.Real.Arithmetic.Elementary++import qualified Data.Number.ER.RnToRm.Approx as FA+import qualified Data.Number.ER.RnToRm.UnitDom.Approx as UFA+import qualified Data.Number.ER.RnToRm.UnitDom.Base as UFB+import qualified Data.Number.ER.Real.Approx as RA+import qualified Data.Number.ER.Real.Approx.Elementary as RAEL++import qualified Data.Number.ER.Real.DomainBox as DBox+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)+import Data.Number.ER.BasicTypes++import Data.Number.ER.Misc++import qualified Data.Map as Map++import Data.Typeable+import Data.Generics.Basics+import Data.Binary++{- only for testing in ghci, to be removed: -}+--import Data.Number.ER.Real.DefaultRepr+--import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom+--import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.PolynomBase+--type FAPU = ERFnInterval (ERChebPoly B) IRA+--fapuConst1 = (UFA.const 0 [1]) :: FAPU+--fapuConst2 = (UFA.const 0 [2]) :: FAPU+{- end of testing specific code -}++data ERFnInterval fb ra =+ ERFnIntervalAny + {+ erfnContext :: ERFnContext+ }+ |+ ERFnInterval + {+ erfnUpper :: fb,+ erfnLowerNeg :: fb,+ erfnContext :: ERFnContext,+ erfnGlobal :: ra+ }+ deriving (Typeable, Data)++instance (Binary a, Binary b) => Binary (ERFnInterval a b) where+ put (ERFnIntervalAny a) = putWord8 0 >> put a+ put (ERFnInterval a b c d) = putWord8 1 >> put a >> put b >> put c >> put d+ get = do+ tag_ <- getWord8+ case tag_ of+ 0 -> get >>= \a -> return (ERFnIntervalAny a)+ 1 -> get >>= \a -> get >>= \b -> get >>= \c -> get >>= \d -> return (ERFnInterval a b c d)+ _ -> fail "no parse"+ ++data ERFnContext =+ ERFnContext+ {+ erfnMaxDegree :: Int,+ erfnCoeffGranularity :: Granularity+ }+ deriving (Show, Typeable, Data)+ +instance Binary ERFnContext where+ put (ERFnContext a b) = put a >> put b+ get = get >>= \a -> get >>= \b -> return (ERFnContext a b)+ + +erfnContextDefault =+ ERFnContext+ {+ erfnMaxDegree = 2,+ erfnCoeffGranularity = 10+ }+ +erfnContextUnify (ERFnContext dg1 gr1) (ERFnContext dg2 gr2) =+ ERFnContext (max dg1 dg2) (max gr1 gr2)++ +instance + (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ Show (ERFnInterval fb ra)+ where+ show (ERFnIntervalAny _) = "ERFnIntervalAny"+ show (ERFnInterval h ln ctxt gl) =+ "\nERFnInterval"+ ++ "\n upper = " ++ show h+ ++ "\n lower = " ++ show (-ln)+-- ++ " global = " ++ show gl ++ "\n"+-- ++ " context = " ++ show ctxt ++ "\n"++instance+ (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ Eq (ERFnInterval fb ra)+ where+ (ERFnInterval h1 ln1 ctxt1 gl1) + == (ERFnInterval h2 ln2 ctxt2 gl2) =+ error "ERFnInterval: equality not implemented yet"+ _ == _ =+ error "ERFnInterval: equality not implemented yet"++instance + (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ Ord (ERFnInterval fb ra) + where+ compare + (ERFnInterval h1 ln1 ctxt1 gl1) + (ERFnInterval h2 ln2 ctxt2 gl2) =+ error "ERFnInterval: comparison not implemented yet"+ compare _ _ =+ error "ERFnInterval: comparison not implemented yet"+ +instance + (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ Num (ERFnInterval fb ra)+ where+ fromInteger n = UFA.const [fromInteger n]+ negate f@(ERFnIntervalAny _) = f+ negate (ERFnInterval h ln ctxt gl) =+ (ERFnInterval ln h ctxt (negate gl))+ (ERFnInterval h1 ln1 ctxt1 gl1) + (ERFnInterval h2 ln2 ctxt2 gl2) =+ ERFnInterval (h1 + h2) (ln1 + ln2) ctxt (gl1 + gl2)+ where+ ctxt = erfnContextUnify ctxt1 ctxt2+ f1 + f2 = ERFnIntervalAny ctxt+ where+ ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)+ (ERFnInterval h1 ln1 ctxt1 gl1) * (ERFnInterval h2 ln2 ctxt2 gl2) =+ ERFnInterval h ln ctxt (gl1 * gl2)+ where+ (h, ln) =+ case (RA.leqReals 0 gl1, RA.leqReals gl1 0, RA.leqReals 0 gl2, RA.leqReals gl2 0) of+ (Just True, _, Just True, _) -> -- both non-negative+ (h1h2, l1l2Neg)+ (_, Just True, _, Just True) -> -- both non-positive+ (l1l2, h1h2Neg)+ (Just True, _, _, Just True) -> -- first non-negative, second non-positive+ (l1h2, h1l2Neg)+ (_, Just True, Just True, _) -> -- first non-positive, second non-negative+ (h1l2, l1h2Neg)+ _ -> -- one of both may be crossing zero+ ((h1h2 `maxP` l1l2) `maxP` (h1l2 `maxP` l1h2),+ (h1h2Neg `maxP` l1l2Neg) `maxP` (h1l2Neg `maxP` l1h2Neg))+ where+ h1h2 = UFB.reduceDegreeUp maxDegr $ h1 * h2+ h1h2Neg = UFB.reduceDegreeUp maxDegr $ (negate h1) * h2+ l1l2 = UFB.reduceDegreeUp maxDegr $ ln1 * ln2+ l1l2Neg = UFB.reduceDegreeUp maxDegr $ (negate ln1) * ln2+ h1l2 = UFB.reduceDegreeUp maxDegr $ h1 * (negate ln2)+ h1l2Neg = UFB.reduceDegreeUp maxDegr $ h1 * ln2+ l1h2 = UFB.reduceDegreeUp maxDegr $ (negate ln1) * h2+ l1h2Neg = UFB.reduceDegreeUp maxDegr $ ln1 * h2+ maxP p1 p2 = fst $ UFB.max maxDegr p1 p2+ + ctxt = erfnContextUnify ctxt1 ctxt2+ maxDegr = erfnMaxDegree ctxt+ f1 * f2 = ERFnIntervalAny ctxt+ where+ ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)+ +instance + (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ Fractional (ERFnInterval fb ra) + where+ fromRational r = UFA.const [fromRational r]+ recip f@(ERFnIntervalAny _) = f+ recip (ERFnInterval h ln ctxt gl) + | certainNoZero =+ ERFnInterval lRecipUp hnRecipUp ctxt (recip gl)+ | otherwise = ERFnIntervalAny ctxt+ where+ certainNoZero =+ certainAboveZero || certainBelowZero+ certainAboveZero =+ UFB.upperBound ix ln < 0+ certainBelowZero = + UFB.upperBound ix h < 0 + hnRecipUp =+ UFB.recipUp maxDegr ix (negate h) + lRecipUp =+ UFB.recipUp maxDegr ix (negate ln)+ maxDegr = erfnMaxDegree ctxt+ ix = int2effIx $ 3 * maxDegr ++instance+ (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ RA.ERApprox (ERFnInterval fb ra) + where+ getGranularity (ERFnIntervalAny ctxt) = erfnCoeffGranularity ctxt+ getGranularity (ERFnInterval h ln ctxt gl) =+ max (erfnCoeffGranularity ctxt) $ + max (UFB.getGranularity h) (UFB.getGranularity ln)+ setGranularity gran (ERFnIntervalAny ctxt) = + ERFnIntervalAny $ ctxt { erfnCoeffGranularity = gran }+ setGranularity gran (ERFnInterval h ln ctxt gl) =+ ERFnInterval + (UFB.setGranularity gran h) (UFB.setGranularity gran ln) + (ctxt { erfnCoeffGranularity = gran }) gl+ setMinGranularity gran (ERFnIntervalAny ctxt) = + ERFnIntervalAny+ (ctxt { erfnCoeffGranularity = max gran (erfnCoeffGranularity ctxt) })+ setMinGranularity gran (ERFnInterval h ln ctxt gl) =+ ERFnInterval + (UFB.setMinGranularity gran h) (UFB.setMinGranularity gran ln) + (ctxt { erfnCoeffGranularity = max gran (erfnCoeffGranularity ctxt) }) gl+-- getPrecision (ERFnIntervalAny _) = 0+-- getPrecision f = intLog 2 (1 + (fst $ RA.integerBounds (FA.volume f))) -- wrong! + (ERFnInterval h1 ln1 ctxt1 gl1) /\ (ERFnInterval h2 ln2 ctxt2 gl2) =+ ERFnInterval (snd $ UFB.min maxDegr h1 h2) (snd $ UFB.min maxDegr ln1 ln2) ctxt (gl1 RA./\ gl2)+ where+ ctxt = erfnContextUnify ctxt1 ctxt2+ maxDegr = erfnMaxDegree ctxt+ (ERFnIntervalAny ctxt1) /\ (ERFnInterval h2 ln2 ctxt2 gl2) =+ ERFnInterval h2 ln2 ctxt gl2+ where+ ctxt = erfnContextUnify ctxt1 ctxt2+ (ERFnInterval h1 ln1 ctxt1 gl1) /\ (ERFnIntervalAny ctxt2) =+ ERFnInterval h1 ln1 ctxt gl1+ where+ ctxt = erfnContextUnify ctxt1 ctxt2+ f1 /\ f2 = ERFnIntervalAny ctxt+ where+ ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)+ leqReals = erfnintLeq++erfnintLeq left right+ | left `isClearlyBelow` right = Just True+ | right `isClearlyStrictlyBelow` left = Just False+ | otherwise = Nothing+ where+ isClearlyBelow (ERFnIntervalAny _) _ = False+ isClearlyBelow _ (ERFnIntervalAny _) = False+ isClearlyBelow f g+ | UFB.upperBound 10 (erfnUpper f + erfnLowerNeg g) <= 0 = True+ | otherwise = False+ isClearlyStrictlyBelow (ERFnIntervalAny _) _ = False+ isClearlyStrictlyBelow _ (ERFnIntervalAny _) = False+ isClearlyStrictlyBelow f g+ | UFB.upperBound 10 (erfnUpper f + erfnLowerNeg g) < 0 = True+ | otherwise = False++instance+ (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ RA.ERIntApprox (ERFnInterval fb ra) + where+-- doubleBounds = :: ira -> (Double, Double) +-- floatBounds :: ira -> (Float, Float)+-- integerBounds :: ira -> (ExtendedInteger, ExtendedInteger)+ bisectDomain maybePt (ERFnIntervalAny c) =+ error "ERFnInterval: RA.bisectDomain: cannot bisect ERFnIntervalAny"+ bisectDomain maybePt (ERFnInterval u ln c g) =+ (ERFnInterval midUp ln c g,+ ERFnInterval u (negate midDown) c g)+ where+ (midDown, midUp) =+ case maybePt of+ Nothing ->+ (negate $ (ln - u) / 2, (u - ln) / 2)+ Just (ERFnInterval uPt lnPt _ _) ->+ (negate lnPt, uPt)+ bounds (ERFnIntervalAny c) =+ error "ERFnInterval: RA.bounds: cannot get bounds for ERFnIntervalAny"+ bounds (ERFnInterval u ln c g) =+ (ERFnInterval (negate ln) ln c g,+ ERFnInterval u (negate u) c g) + (ERFnInterval u1 ln1 c1 g1) \/ (ERFnInterval u2 ln2 c2 g2) =+ ERFnInterval u ln c (g1 RA.\/ g2)+ where+ u = UFB.maxUp maxDegree u1 u2+ ln = UFB.maxUp maxDegree ln1 ln2+ c = erfnContextUnify c1 c2+ maxDegree = erfnMaxDegree c+ (ERFnIntervalAny ctxt1) \/ (ERFnInterval h2 ln2 ctxt2 gl2) =+ ERFnIntervalAny ctxt+ where+ ctxt = erfnContextUnify ctxt1 ctxt2+ (ERFnInterval h1 ln1 ctxt1 gl1) \/ (ERFnIntervalAny ctxt2) =+ ERFnIntervalAny ctxt+ where+ ctxt = erfnContextUnify ctxt1 ctxt2+ f1 \/ f2 = ERFnIntervalAny ctxt+ where+ ctxt = erfnContextUnify (erfnContext f1) (erfnContext f2)++instance+ (UFB.ERUnitFnBase boxb boxra varid b ra fb, RAEL.ERApproxElementary ra) =>+ RAEL.ERApproxElementary (ERFnInterval fb ra) + where+ -- default abs does not work because we do not have Prelude.abs+ abs _ f@(ERFnIntervalAny _) = f+ abs _ (ERFnInterval u ln c g) =+ ERFnInterval maxulnUp maxunl0Dn c (abs g)+ where+ maxDegree = erfnMaxDegree c+ maxulnUp = snd $ UFB.max maxDegree u ln + maxunl0Dn =+ fst $ UFB.max maxDegree 0 $+ fst $ UFB.max maxDegree (- u) (- ln)+ exp ix f@(ERFnIntervalAny _) = f+ exp ix (ERFnInterval u ln c g) =+ ERFnInterval uExp lExpNeg c (RAEL.exp ix g)+ where+ maxDegree = erfnMaxDegree c+-- ix = int2effIx maxDegree+ uExp = snd $ UFB.exp maxDegree ix u+ lExpNeg = + negate $ fst $ UFB.exp maxDegree ix (negate ln) + sin ix f@(ERFnIntervalAny c) = + ERFnInterval 1 1 c ((-1) RA.\/ 1)+ sin ix (ERFnInterval u ln c g) =+-- unsafePrint+-- (+-- "ERFnInterval: RAEL.sin: "+-- ++ "\n u = " ++ show u+-- ++ "\n ln = " ++ show ln+-- ++ "\n uSin = " ++ show uSin+-- ++ "\n lSinNeg = " ++ show lSinNeg+-- ) $+ ERFnInterval uSin lSinNeg c (RAEL.sin ix g)+ where+ maxDegree = erfnMaxDegree c+-- ix = int2effIx maxDegree+ uSin = snd $ UFB.sin maxDegree ix u+ lSinNeg = + negate $ fst $ UFB.sin maxDegree ix (negate ln) + cos ix f@(ERFnIntervalAny c) =+ ERFnInterval 1 1 c ((-1) RA.\/ 1)+ cos ix (ERFnInterval u ln c g) =+-- unsafePrint+-- (+-- "ERFnInterval: RAEL.sin: "+-- ++ "\n u = " ++ show u+-- ++ "\n ln = " ++ show ln+-- ++ "\n uSin = " ++ show uSin+-- ++ "\n lSinNeg = " ++ show lSinNeg+-- ) $+ ERFnInterval uCos lCosNeg c (RAEL.cos ix g)+ where+ maxDegree = erfnMaxDegree c+-- ix = int2effIx maxDegree+ uCos = snd $ UFB.cos maxDegree ix u+ lCosNeg = + negate $ fst $ UFB.cos maxDegree ix (negate ln) ++instance + (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ FA.ERFnApprox boxra varid ra ra (ERFnInterval fb ra)+ where+ check callerLocation f@(ERFnIntervalAny c) = f+ check callerLocation (ERFnInterval u ln c g) =+ ERFnInterval + (UFB.check (callerLocation ++ "upper: ") u) + (UFB.check (callerLocation ++ "neg lower: ") ln) + c g + domra2ranra _ = id+ ranra2domra _ = id+ setMaxDegree maxDegr (ERFnIntervalAny c) =+ ERFnIntervalAny (c { erfnMaxDegree = maxDegr } )+ setMaxDegree maxDegr (ERFnInterval u ln c g) =+ ERFnInterval + (UFB.reduceDegreeUp maxDegr u)+ (UFB.reduceDegreeUp maxDegr ln)+ (c { erfnMaxDegree = maxDegr } )+ g+ getRangeApprox (ERFnIntervalAny _) = RA.bottomApprox + getRangeApprox (ERFnInterval u ln c g) =+ UFB.raFromEndpoints u+ (+ (UFB.upperBound 10 u)+ ,+ (- (UFB.upperBound 10 ln))+ )+ scale ratio f@(ERFnIntervalAny c) = f+ scale ratio f@(ERFnInterval u ln c g) = + case RA.compareReals ratio 0 of+ Just GT -> + ERFnInterval (UFB.scaleApproxUp ratio u) (UFB.scaleApproxUp ratio ln) c g+ Just LT -> + ERFnInterval (UFB.scaleApproxUp (- ratio) ln) (UFB.scaleApproxUp (- ratio) u) c g+ _ -> + (UFA.const [ratio]) * f+ eval ptBox (ERFnIntervalAny c) = [RA.bottomApprox]+ eval ptBox (ERFnInterval u ln c g) =+ [lo RA.\/ up]+ where+ up = UFB.evalApprox ptBox u+ lo = negate $ UFB.evalApprox ptBox ln+ partialEval substitutions f@(ERFnIntervalAny c) = f+ partialEval substitutions (ERFnInterval u ln c g) =+ (ERFnInterval uP lnP c g)+ where+ uP = UFB.partialEvalApproxUp substitutions u+ lnP = UFB.partialEvalApproxUp substitutions ln++ composeThin+ f@(ERFnIntervalAny ctxt)+ substitutions =+ f+ composeThin+ f@(ERFnInterval h1 ln1 ctxt1 gl1)+ substitutions =+ (ERFnInterval h ln ctxt1 gl1)+ where+ h = UFB.composeUp maxDegree h1 ufbSubstitutions + ln = UFB.composeUp maxDegree ln1 ufbSubstitutions+ ufbSubstitutions = Map.map erfnUpper substitutions+ maxDegree = erfnMaxDegree ctxt1 +-- ctxt = erfnContextUnify ctxt1 ctxt2++instance + (UFB.ERUnitFnBase boxb boxra varid b ra fb) =>+ UFA.ERUnitFnApprox boxra varid ra ra (ERFnInterval fb ra)+ where+ bottomApprox =+ ERFnIntervalAny erfnContextDefault+ const [val] + | RA.isBounded val =+ ERFnInterval+ {+ erfnUpper = fbH,+ erfnLowerNeg = fbLNeg,+ erfnContext = context,+ erfnGlobal = val+ }+ | otherwise =+ ERFnIntervalAny context + where+ fbH = UFB.const valH+ fbLNeg = UFB.const (negate valL)+ (valL, valH) = UFB.raEndpoints fbH val+ context = + erfnContextDefault+ {+ erfnCoeffGranularity = RA.getGranularity val+ }+ affine [val] coeffsSingletons+ | RA.isBounded val && (and $ map (RA.isBounded . head) $ Map.elems coeffsSingletons) =+ ERFnInterval+ {+ erfnUpper = fbH,+ erfnLowerNeg = fbLNeg,+ erfnContext = context,+ erfnGlobal = + UFB.raFromEndpoints fbH+ (valL - coeffCorr - coeffsAbsSum, + valH + coeffCorr + coeffsAbsSum)+ }+ | otherwise =+ ERFnIntervalAny context+ where+ coeffs = Map.map (\[a] -> a) coeffsSingletons+ coeffGranularity =+ Map.fold max (RA.getGranularity val) (Map.map RA.getGranularity coeffs)+ coeffsMsCorrs = + Map.map (\(l,h) ->+ (B.setMinGranularity coeffGranularity (l + h)/2, + B.setMinGranularity coeffGranularity (h - l)/2)) $+ Map.map (UFB.raEndpoints fbH) $ coeffs+ coeffCorr = Map.fold (+) 0 $ Map.map snd coeffsMsCorrs+ coeffsAbsSum = Map.fold (+) 0 $ Map.map (abs . fst) coeffsMsCorrs+ fbH = UFB.affine (valH + coeffCorr) (Map.map fst coeffsMsCorrs)+ fbLNeg = UFB.affine (negate (valL - coeffCorr)) (Map.map (negate . fst) coeffsMsCorrs)+ (valL, valH) = UFB.raEndpoints fbH val+ context = + erfnContextDefault+ {+ erfnCoeffGranularity = coeffGranularity+ }+ intersectMeasureImprovement ix vars+ f1@(ERFnIntervalAny ctxt1) + f2@(ERFnIntervalAny ctxt2) =+ (ERFnIntervalAny ctxt, RA.bottomApprox)+ where+ ctxt = erfnContextUnify ctxt1 ctxt2+ intersectMeasureImprovement ix vars+ f1@(ERFnIntervalAny ctxt1) + f2@(ERFnInterval h2 ln2 ctxt2 gl2) =+ (ERFnInterval h2 ln2 ctxt gl2, 1 / 0)+ where+ ctxt = erfnContextUnify ctxt1 ctxt2+ intersectMeasureImprovement ix vars+ f1@(ERFnInterval h1 ln1 ctxt1 gl1) + f2@(ERFnIntervalAny ctxt2) = + (ERFnInterval h1 ln1 ctxt gl1, 1)+ where+ ctxt = erfnContextUnify ctxt1 ctxt2+ intersectMeasureImprovement ix vars+ f1@(ERFnInterval h1 ln1 ctxt1 gl1) + f2@(ERFnInterval h2 ln2 ctxt2 gl2) =+ case RA.compareReals improvementRA 1 of+ Just LT -> (f1, 1) -- intersection made it worse, keep original+ _ -> (intersection, improvementRA)+ where+ intersection = f1 RA./\ f2+ improvementRA + | 0 `RA.refines` intersectionVolume && 0 `RA.refines` f1Volume = 1+-- error $ +-- "ERFnInterval: intersectMeasureImprovement: inconsistent result: " +-- ++ show intersection+ | otherwise =+ f1Volume / intersectionVolume+ intersectionVolume = UFA.volume vars intersection+ f1Volume = UFA.volume vars f1+ ctxt = erfnContextUnify ctxt1 ctxt2+ volume vars (ERFnIntervalAny c) = 1/0+ volume vars (ERFnInterval u ln c g) =+-- unsafePrint ("ERFnInterval: volume: result = " ++ show result) $ result+-- where+-- result =+ UFB.raFromEndpoints u $ UFB.volumeAboveZero vars (u + ln)+ integrate + ix (ERFnInterval u ln c g) x + origin (ERFnInterval uInit lnInit cInit gInit) =+-- unsafePrint+-- (+-- "ERFnInterval: integrate: " +-- ++ "\n u = " ++ show u+-- ++ "\n ln = " ++ show ln+-- ++ "\n origin = " ++ show origin+-- ++ "\n uInit = " ++ show uInit+-- ++ "\n lnInit = " ++ show lnInit+-- ++ "\n uIuL = " ++ show uIuL+-- ++ "\n uIuU = " ++ show uIuU+-- ++ "\n uIuOriginL = " ++ show uIuOriginL+-- ++ "\n uIuOriginU = " ++ show uIuOriginU+-- ++ "\n lnIuL = " ++ show lnIuL+-- ++ "\n lnIuU = " ++ show lnIuU+-- ++ "\n lnIuOriginL = " ++ show lnIuOriginL+-- ++ "\n lnIuOriginU = " ++ show lnIuOriginU+-- ++ "\n uIov = " ++ show uIov+-- ++ "\n lnIov = " ++ show lnIov+-- )+ (ERFnInterval uIov lnIov c gIov)+ where+ -- perform raw integration of both bounds:+ (uIuL, uIuU) = +-- mapPair (UFB.reduceDegreeDown maxDegree, UFB.reduceDegreeUp maxDegree) $ + UFB.integrate x u+ (lnIuL, lnIuU) = +-- mapPair (UFB.reduceDegreeDown maxDegree, UFB.reduceDegreeUp maxDegree) $ + UFB.integrate x ln+ maxDegree = erfnMaxDegree c+ -- constrain the raw integrals to the origin:+ uIuOriginL = UFB.composeDown maxDegree uIuL substXOrigin+ uIuOriginU = UFB.composeUp maxDegree uIuU substXOrigin+ lnIuOriginL = UFB.composeDown maxDegree lnIuL substXOrigin+ lnIuOriginU = UFB.composeUp maxDegree lnIuU substXOrigin+ substXOrigin = Map.singleton x originUFB+ originUFB = UFB.const $ fst $ UFB.raEndpoints u origin+ -- adjust the raw integrated functions enclose the initial condition function: + uIov = + UFB.reduceDegreeUp maxDegree $+ uIuU + uInit - uIuOriginL + (uIuOriginU - uIuOriginL)+ lnIov = + UFB.reduceDegreeUp maxDegree $+ lnIuU + lnInit - lnIuOriginL + (lnIuOriginU - lnIuOriginL)+ + gIov = + gInit + g * ((1 - origin) RA.\/ (-1 - origin))+
+ src/Data/Number/ER/RnToRm/UnitDom/Base.hs view
@@ -0,0 +1,348 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.Base+ Description : class abstracting imprecise function arithmetic on [-1,1]^n+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++ A class abstracting function arithmetic with directed rounding.+ It is used to describe a boundary for an approximation+ to a real function on the interval [-1,1]^n.+ + To be imported qualified, usually with the synonym UFB.+-}+module Data.Number.ER.RnToRm.UnitDom.Base where++import Prelude hiding (min, max, recip)++import qualified Data.Number.ER.Real.DomainBox as DBox+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)+import Data.Number.ER.BasicTypes+import qualified Data.Number.ER.Real.Base as B+import qualified Data.Number.ER.Real.Approx as RA++import qualified Data.Map as Map++import Data.Typeable++class + (B.ERRealBase b, RA.ERIntApprox ra, Fractional ufb, Ord ufb,+ DomainBox boxb varid b, DomainIntBox boxra varid ra) => + ERUnitFnBase boxb boxra varid b ra ufb+ | ufb -> boxb boxra varid b ra+ where+ {-| + Check internal consistency of the function and report problem if any.+ -}+ check :: + String {-^ indentification of caller location for easier debugging -} -> + ufb -> ufb+ getGranularity :: ufb -> Granularity+ setMinGranularity :: Granularity -> ufb -> ufb+ setGranularity :: Granularity -> ufb -> ufb+ {-| Construct a constant function. -}+ const :: b -> ufb+ {-| Construct an affine function. -}+ affine :: + b {-^ value at 0 -} ->+ Map.Map varid b {-^ ascent of each base vector -} -> + ufb+ {-| + Multiply a function by a scalar, + rounding downwards and upwards. + -} + scale :: b -> ufb -> (ufb, ufb) + {-| + Multiply a function by an approximation of a scalar, + rounding downwards and upwards. + -} + scaleApprox :: ra -> ufb -> (ufb, ufb) + {-| + Multiply a function by an approximation of a scalar, + rounding downwards. + -} + scaleApproxDown :: ra -> ufb -> ufb+ scaleApproxDown ratio = fst . scaleApprox ratio + {-| + Multiply a function by an approximation of a scalar, + rounding upwards. + -} + scaleApproxUp :: ra -> ufb -> ufb+ scaleApproxUp ratio = snd . scaleApprox ratio + {-| + Get the degree of this particular function.+ + If the function is a polynomial, this function should+ return its degree. + -}+ getDegree :: ufb -> Int+ {-| + Decrease the degree of function approximation, + rounding pointwise downwards and upwards.+ -}+ reduceDegree :: Int -> ufb -> (ufb, ufb)+ {-| + Decrease the degree of function approximation, rounding pointwise downwards.+ -}+ reduceDegreeDown :: Int -> ufb -> ufb+ reduceDegreeDown maxDegr = fst . reduceDegree maxDegr+ {-| + Decrease the degree of function approximation, rounding pointwise upwards.+ -}+ reduceDegreeUp :: Int -> ufb -> ufb+ reduceDegreeUp maxDegr = snd . reduceDegree maxDegr+ {-| + Approximate the integral of p (with 0 at 0) from below and from above.+ -}+ integrate :: + varid {-^ variable to integrate by -} -> + ufb {-^ p(x) -} -> + (ufb, ufb)+ {-| Approximate the integral of p (with 0 at 0) from below. -}+ integrateDown :: + varid {-^ variable to integrate by -} -> + ufb {-^ p(x) -} -> + ufb+ integrateDown x = fst . integrate x+ {-| Approximate the integral of p (with 0 at 0) from above. -}+ integrateUp :: + varid {-^ variable to integrate by -} -> + ufb {-^ p(x) -} -> + ufb+ integrateUp x = snd . integrate x+ {-| + Measure the volume between a function + and the zero hyperplane on the domain @[-1,1]^n@.+ -}+ volumeAboveZero :: + [varid] {-^ axes to include in the measuring domain -} -> + ufb -> (b,b)+ {-|+ Find an upper bound of the function over @[-1,1]^n@.+ -}+ upperBound :: EffortIndex -> ufb -> b+ {-|+ Find a lower bound of the function over @[-1,1]^n@.+ -}+ lowerBound :: EffortIndex -> ufb -> b+ lowerBound ix f = negate $ upperBound ix (negate f)+ {-| + Approximate the function max(0,p(x)) from below and from above.+ -}+ nonneg ::+ Int {-^ max degree for result -} -> + ufb {-^ p(x) -} -> + (ufb, ufb)+ {-| + Approximate the function 1/p(x) from below and from above.+ -}+ recip :: + Int {-^ max degree for result -} ->+ EffortIndex -> + ufb {-^ p(x) -} -> + (ufb, ufb)+ {-| + Approximate the function 1/p(x) from below.+ -}+ recipDown :: Int -> EffortIndex -> ufb -> ufb+ recipDown maxDegr ix a = fst $ recip maxDegr ix a+ {-| + Approximate the function 1/p(x) from above.+ -}+ recipUp :: Int -> EffortIndex -> ufb -> ufb+ recipUp maxDegr ix a = snd $ recip maxDegr ix a+ {-|+ Approximate the function max(p_1(x),p_2(x)) from below and from above.+ -}+ max :: + Int {-^ max degree for result -} -> + ufb {-^ p_1(x) -} -> + ufb {-^ p_2(x) -} -> + (ufb, ufb)+ {-|+ Approximate the function max(p_1(x),p_2(x)) from below.+ -}+ maxDown :: + Int {-^ max degree for result -} -> + ufb {-^ p_1(x) -} -> + ufb {-^ p_2(x) -} -> + ufb+ maxDown maxDegr a b = fst $ max maxDegr a b+ {-|+ Approximate the function max(p_1(x),p_2(x)) from above.+ -}+ maxUp :: + Int {-^ max degree for result -} -> + ufb {-^ p_1(x) -} -> + ufb {-^ p_2(x) -} -> + ufb+ maxUp maxDegr a b = snd $ max maxDegr a b+ {-|+ Approximate the function min(p_1(x),p_2(x)) from below and from above.+ -}+ min :: + Int {-^ max degree for result -} -> + ufb {-^ p_1(x) -} -> + ufb {-^ p_2(x) -} -> + (ufb, ufb)+ min maxDegr p1 p2 = -- default implementation using symmetry with ufbMax+ (negate hi, negate lo)+ where+ (lo, hi) = max maxDegr (negate p1) (negate p2)+ {-|+ Approximate the function min(p_1(x),p_2(x)) from below.+ -}+ minDown :: + Int {-^ max degree for result -} -> + ufb {-^ p_1(x) -} -> + ufb {-^ p_2(x) -} -> + ufb+ minDown maxDegr a b = fst $ min maxDegr a b+ {-|+ Approximate the function min(p_1(x),p_2(x)) from above.+ -}+ minUp :: + Int {-^ max degree for result -} -> + ufb {-^ p_1(x) -} -> + ufb {-^ p_2(x) -} -> + ufb+ minUp maxDegr a b = snd $ min maxDegr a b+ {-|+ Approximate @sqrt(p(x))@ from below and from above.+ -}+ sqrt :: + Int {-^ max degree for result -} -> + EffortIndex {-^ how hard to try when approximating exp as a polynomial -} -> + ufb {-^ p(x) -} -> + (ufb, ufb)+ {-|+ Approximate @exp(p(x))@ from below and from above.+ -}+ exp :: + Int {-^ max degree for result -} -> + EffortIndex {-^ how hard to try when approximating exp as a polynomial -} -> + ufb {-^ p(x) -} -> + (ufb, ufb)+ {-| + Approximate @log(p(x))@ from below and from above.+ -}+ log :: + Int {-^ max degree for result -} -> + EffortIndex {-^ how hard to try when approximating log as a polynomial -} -> + ufb {-^ p(x) -} -> + (ufb, ufb)+ {-| + Approximate @sin(p(x))@ from below and from above.+ -}+ sin :: + Int {-^ max degree for result -} -> + EffortIndex {-^ how hard to try when approximating sin as a polynomial -} -> + ufb {-^ p(x) -} -> + (ufb, ufb)+ {-|+ Approximate @cos(p(x))@ from below and from above.+ -}+ cos :: + Int {-^ max degree for result -} -> + EffortIndex {-^ how hard to try when approximating cos as a polynomial -} -> + ufb {-^ p(x) -} -> + (ufb, ufb)+ {-| + Evaluate at a point, rounding upwards and downwards.+ -}+ eval :: boxb -> ufb -> (b, b)+ {-| + Evaluate at a point, rounding downwards.+ -}+ evalDown :: boxb -> ufb -> b+ evalDown pt = fst . eval pt+ {-| + Evaluate at a point, rounding downwards.+ -}+ evalUp :: boxb -> ufb -> b+ evalUp pt = snd . eval pt+ {-|+ Safely evaluate at a point using a real number approximation+ for both the point and the result.+ -}+ evalApprox :: boxra -> ufb -> ra+ {-|+ Partially evaluate at a lower-dimensional point + given using a real number approximation.+ Approximate the resulting function from below and from above.+ -}+ partialEvalApprox :: boxra -> ufb -> (ufb, ufb)+ {-|+ Partially evaluate at a lower-dimensional point + given using a real number approximation.+ Approximate the resulting function from below.+ -}+ partialEvalApproxDown :: boxra -> ufb -> ufb+ partialEvalApproxDown substitutions = fst . partialEvalApprox substitutions+ {-|+ Partially evaluate at a lower-dimensional point + given using a real number approximation.+ Approximate the resulting function from above.+ -}+ partialEvalApproxUp :: boxra -> ufb -> ufb+ partialEvalApproxUp substitutions = snd . partialEvalApprox substitutions+ {-| + Compose two functions, rounding upwards and downwards+ provided each @f_v@ ranges within the domain @[-1,1]@. + -} + compose ::+ Int {-^ max degree for result -} -> + ufb {-^ function @f@ -} -> + Map.Map varid ufb + {-^ variables to substitute and for each variable @v@, + function @f_v@ to substitute for @v@ + that maps @[-1,1]@ into @[-1,1]@ -} ->+ (ufb, ufb) {-^ upper and lower bounds of @f[v |-> f_v]@ -}+ {-| + Compose two functions, rounding downwards+ provided each @f_v@ ranges within the domain @[-1,1]@. + -} + composeDown ::+ Int {-^ max degree for result -} -> + ufb {-^ function @f1@ -} -> + Map.Map varid ufb + {-^ variables to substitute and for each variable @v@, + function @f_v@ to substitute for @v@ + that maps @[-1,1]@ into @[-1,1]@ -} ->+ ufb {-^ a lower bound of @f1.f2@ -}+ composeDown maxDegr f = fst . compose maxDegr f + {-| + Compose two functions, rounding upwards+ provided each @f_v@ ranges within the domain @[-1,1]@. + -} + composeUp ::+ Int {-^ max degree for result -} -> + ufb {-^ function @f1@ -} -> + Map.Map varid ufb + {-^ variables to substitute and for each variable @v@, + function @f_v@ to substitute for @v@ + that maps @[-1,1]@ into @[-1,1]@ -} ->+ ufb {-^ an upper bound of @f1.f2@ -}+ composeUp maxDegr f = snd . compose maxDegr f + {-|+ Convert from the interval type to the base type.+ (The types are determined by the given example function.)+ -}+ raEndpoints :: + ufb {-^ this parameter is not used except for type checking -} -> + ra -> + (b,b)+ {-|+ Convert from the base type to the interval type. + (The types are determined by the given example function.)+ -}+ raFromEndpoints :: + ufb {-^ this parameter is not used except for type checking -} -> + (b,b) ->+ ra+
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs view
@@ -0,0 +1,93 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE UndecidableInstances #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom+ Description : polynoms in the Chebyshev basis of the 1st kind+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Arithmetic of multivariate polynomials + represented by their coefficients it the Chebyshev basis.+ + The polynomials are never to be used outside the domain @[-1,1]^n@.+ + All operations are rounded in such a way that the resulting polynomial+ is a /point-wise upper or lower bound/ of the exact result. +-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom+(+ ERChebPoly(..), TermKey+) +where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field+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.Integration+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary++import qualified Data.Number.ER.RnToRm.UnitDom.Base as UFB+import qualified Data.Number.ER.Real.Base as B+import Data.Number.ER.Real.Approx.Interval+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)++{- code for testing purpose, to be deleted later -}+import Data.Number.ER.Real.DefaultRepr+import Data.Number.ER.Real.DomainBox.IntMap+type P = ERChebPoly (Box Int) B+x0 = chplVar 0 :: P+x1 = chplVar 1 :: P+x2 = chplVar 2 :: P+x3 = chplVar 3 :: P+x4 = chplVar 4 :: P+p1 = x1 * x1 * x1 + x1 * (x2 + 2) * (x3 - 3)+{- end of code for testing purposes -}+++instance + (B.ERRealBase rb, RealFrac rb,+ DomainBox box varid Int, Ord box,+ DomainBoxMappable boxb boxbb varid rb [(rb,rb)],+ DomainBoxMappable boxra boxras varid (ERInterval rb) [ERInterval rb],+ DomainIntBox boxra varid (ERInterval rb)) =>+ (UFB.ERUnitFnBase boxb boxra varid rb (ERInterval rb) (ERChebPoly box rb))+ where+ check = chplCheck+ getGranularity = chplGetGranularity+ setMinGranularity = chplSetMinGranularity+ setGranularity = chplSetGranularity+ const = chplConst+ affine = chplAffine+ scale = chplScale+ scaleApprox (ERInterval ratioDown ratioUp) = chplScaleApprox (ratioDown, ratioUp) +-- Arity = chplGetArity+ getDegree = chplGetDegree+ reduceDegree = chplReduceDegree+ volumeAboveZero = chplVolumeAboveZero+ integrate = chplIntegrate+ upperBound = chplUpperBoundAffine+-- upperBound = chplUpperBoundQuadr+ nonneg = chplNonneg+ recip = chplRecip+ max = chplMax+ sqrt = chplSqrt+ exp = chplExp+ log = chplLog+ sin = chplSineCosine True+ cos = chplSineCosine False+ eval = chplEval+ evalApprox ufb x = chplEvalApprox (\ b -> ERInterval b b) ufb x+ partialEvalApprox substitutions ufb = + chplPartialEvalApprox (UFB.raEndpoints ufb) substitutions ufb+ raEndpoints _ (ERInterval l h) = (l,h)+ raEndpoints _ ERIntervalAny = (- B.plusInfinity, B.plusInfinity)+ raFromEndpoints _ (l,h) = normaliseERInterval (ERInterval l h)+ compose = chplCompose+
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs view
@@ -0,0 +1,279 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+ Description : (internal) polynomial datatype and simple functions+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Internal module for "Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom".+ + Definition of the polynomial datatype and simple related functions.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic +where++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, DomainIntBox)+import Data.Number.ER.Misc++import qualified Data.Map as Map++import Data.Typeable+import Data.Generics.Basics+import Data.Binary++errorModule msg = error $ "Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom: " ++ msg++{-|+ A polynomial represented by its coefficients it the Chebyshev basis.+ + The polynomials are never to be used outside the domain @[-1,1]^n@.+ + All operations are rounded in such a way that the resulting polynomial+ is a /point-wise upper or lower bound/ of the exact result. +-}+data ERChebPoly box b =+ ERChebPoly+-- Map (MultiSet Int) b+ {+ chplCoeffs :: (Map.Map (TermKey box) b)+ }+ deriving (Eq, Typeable, Data)+ +instance (Ord a, Binary a, Binary b) => Binary (ERChebPoly a b) where+ put (ERChebPoly a) = put a+ get = get >>= \a -> return (ERChebPoly a)+ + +chplCheck prgLocation p@(ERChebPoly coeffs)+ | ok = p+ | otherwise = + unsafePrint (prgLocation ++ " problem with p = \n" ++ show p) p+ where+ ok = + Map.fold (&&) True $ Map.map coeffOK coeffs+ coeffOK c =+ not $ B.isERNaN c+ +type TermKey box = box+ +chplConstTermKey :: (DomainBox box varid d) => box+chplConstTermKey = DBox.noinfo++chplIsConstTermKey :: (DomainBox box varid d) => box -> Bool+chplIsConstTermKey = DBox.isNoinfo++chplTermOrder :: (DomainBox box varid d, Num d) => box -> d+chplTermOrder termKey = DBox.fold (+) 0 termKey++chplTermArity :: (DomainBox box varid d) => box -> Int+chplTermArity termKey = length $ DBox.keys termKey++{-|+ Inspect all terms of the polynomial and return the + degree of the highest degree term.+-}+chplGetDegree ::+ (B.ERRealBase b, DomainBox box varid d, Num d, Ord d) =>+ (ERChebPoly box b) ->+ d+chplGetDegree (ERChebPoly coeffs) =+ foldl max 0 $ map chplTermOrder $ Map.keys coeffs++ +-- chplGetArity = length . chplGetVars + +chplGetVars (ERChebPoly coeffs) =+ DBox.keys $ foldl DBox.union DBox.noinfo $ Map.keys coeffs++chplGetGranularity (ERChebPoly coeffs) =+ foldl max 0 $ map B.getGranularity $ Map.elems coeffs+ +chplSetMinGranularity gran (ERChebPoly coeffs) =+ ERChebPoly $ Map.map (B.setMinGranularity gran) coeffs+ +chplSetGranularity gran (ERChebPoly coeffs) =+ ERChebPoly $ Map.map (B.setGranularity gran) coeffs+ +chplConst :: + (B.ERRealBase b, DomainBox box varid Int, Ord box) => + b -> + ERChebPoly box b+chplConst val = + (ERChebPoly $ Map.singleton chplConstTermKey val)+ +{-|+ make a basic "x" polynomial for a given variable number +-}+chplVar :: + (B.ERRealBase b, DomainBox box varid Int, Ord box) => + varid -> + ERChebPoly box b+chplVar varName =+ ERChebPoly $ Map.singleton (DBox.singleton varName 1) 1++--{-|+-- Make a univariate polynomial given by a series of coefficients+-- in the Chebyshev basis. +---}+--chplMakeUnivariate ::+-- (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+-- varid ->+-- [(Int, b)] {-^ list of pairs: degree of Chebyshev polynomial + coefficient -} ->+-- ERChebPoly box b+--chplMakeUnivariate varName powCoeffPairs =+-- ERChebPoly $ Map.fromList $ map encodePow powCoeffPairs+-- where+-- encodePow (pow, coeff) =+-- (DBox.singleton varName pow, coeff)++chplNormaliseDown, chplNormaliseUp ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ ERChebPoly box b -> ERChebPoly box b+chplNormaliseUp = snd . chplNormalise+chplNormaliseDown = fst . chplNormalise++chplNormalise ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ ERChebPoly box b -> (ERChebPoly box b, ERChebPoly box b)+chplNormalise (ERChebPoly coeffs) =+ (ERChebPoly $ coeffsNo0T0Down,+ ERChebPoly $ coeffsNo0T0Up)+ where+ coeffsNo0T0Down =+ Map.insertWith plusDown chplConstTermKey err coeffsNo0T0+ coeffsNo0T0Up =+ Map.insertWith plusUp chplConstTermKey err coeffsNo0T0+ (coeffsNo0T0, err) = + foldl addTermNo0T0 (Map.empty, 0) $ Map.toList coeffs+ addTermNo0T0 (prevCoeffs, prevErr) (term, coeff) + | coeff == 0 =+ (prevCoeffs, prevErr)+ | otherwise =+ (newCoeffs, newErr)+ where+ newTerm =+ DBox.filter (> 0) term+ newCoeffs = + Map.insert newTerm newCoeffUp prevCoeffs+ newCoeffUp = prevCoeff + coeff+ newCoeffDown = prevCoeff `plusDown` coeff+ prevCoeff =+ Map.findWithDefault 0 newTerm prevCoeffs+ newErr = newCoeffUp - newCoeffDown++instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Show (ERChebPoly box b)+ where+-- show = chplShow True+ show = chplShow False++{-|+ Convert a polynomial to a string representation,+ using the ordinary x^n basis.+-}+chplShow :: + (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ Bool {-^ show the polynomial also in its native Chebyshev basis -} ->+ ERChebPoly box b ->+ String+chplShow showChebyshevBasis (ERChebPoly coeffs) + | showChebyshevBasis = "\n" ++ inChebBasis ++ " = \n" ++ inXBasis+ | otherwise = inXBasis+ where+ inChebBasis = + showCoeffs showTermT $ coeffs+ inXBasis = + showCoeffs showTermX $ chebToXBasis coeffs+ showCoeffs showTerm coeffs =+ concatWith " + " $ map showTerm $ Map.toAscList coeffs+ showTermT (term, coeff)+ | chplIsConstTermKey term = show coeff+ | otherwise = + show coeff ++ "*" ++ (concat $ map showT $ DBox.toList term) + showTermX (term, coeff)+ | chplIsConstTermKey term = showC coeff+ | otherwise = + showC coeff ++ "*" ++ (concat $ map showX $ DBox.toList term) + showT (var, deg) = "T" ++ show deg ++ "(" ++ showVar var ++ ")"+ showX (var, deg) = showVar var ++ "^" ++ show deg+ showC = B.showDiGrCmp 8 False False++{-|+ conversion of polynomials from Chebyshev basis to the X^n basis+ + (not exact - suffering from rounding in the coefficient conversions)+-}+chebToXBasis ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ (Map.Map (TermKey box) b) {-^ polynomial in Chebyshev basis -} ->+ (Map.Map (TermKey box) b) {-^ approxition of the equivalent polynomial in X^n basis -}+chebToXBasis coeffs =+ Map.foldWithKey addTerm Map.empty coeffs+ where+ addTerm term coeff prevXCoeffs =+ Map.unionWith (+) prevXCoeffs $+ Map.map (\ c -> coeff * (fromInteger c)) $ + termXterms term++{-|+ conversion of one Chebyshev term to the X^n basis+-}+termXterms ::+ (DomainBox box varid Int, Ord box) =>+ TermKey box + {-^ a Chebyshev term represented by the Chebyshev degrees + for each variable in the term -} ->+ Map.Map (TermKey box) Integer+ {-^ the polynomial equivalent to the given Chebyshev term + (using integer coefficients) -}+termXterms term =+ foldl addCombination Map.empty $ + allCombinations $ + map (mapSnd $ \ deg -> chebyXCoeffsLists !! deg) $ + DBox.toList term+ where+ addCombination prevMap (varPowerCoeffTriples) =+ Map.insertWith (+) term coeffProduct prevMap+ where+ term = + DBox.fromList $+ filter (\(v,p) -> p > 0) $+ map (\(v,(p,_)) -> (v,p)) varPowerCoeffTriples + coeffProduct =+ fromInteger $+ product $ + map (\(_,(_,c)) -> c) varPowerCoeffTriples+ +{-| Chebyshev polynomials expressed as associative lists power -> coeff -}+chebyXCoeffsLists ::+ (Num d1, Enum d1, Num d2, Enum d2) =>+ [[(d1, d2)]]+chebyXCoeffsLists =+ map convertCoeffs chebyXCoeffs+ where+ convertCoeffs coeffs =+ filter ((/= 0) . snd) $ zip [0,1..] coeffs++{-| Chebyshev polynomials expressed as lists of integer coefficients for powers 0,1,2... -}+chebyXCoeffs ::+ (Num d, Enum d) =>+ [[d]]+chebyXCoeffs =+ aux + [1] -- T_0(x) = 1+ [0,1] -- T_1(x) = x+ where+ aux tnM2 tnM1 =+ tnM2 : (aux tnM1 (newTerm tnM2 tnM1))+ newTerm tnM2 tnM1 =+ zipWith (-) (0 : (map (*2) tnM1)) (tnM2 ++ [0,0..])+ -- T_n(x) = 2 * x * T_{n-1}(x) - T_{n-2}(x)+
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs view
@@ -0,0 +1,293 @@+{-# LANGUAGE FlexibleContexts #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds+ Description : (internal) bounds of single and multiple polynomials+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++ Internal module for "Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom".+ + Implementation of various functions related to the bounds of polynomials. +-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds +where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field++import qualified Data.Number.ER.Real.Approx as RA+import qualified Data.Number.ER.Real.Base as B+import Data.Number.ER.Real.Approx.Interval+import Data.Number.ER.Real.Arithmetic.LinearSolver+import qualified Data.Number.ER.Real.DomainBox as DBox+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox)+import Data.Number.ER.BasicTypes+import Data.Number.ER.Misc++import qualified Data.Map as Map++import Data.List++{-|+ Find an upper bound on a polynomial over the + unit domain [-1,1]^n. +-}+chplUpperBoundAffine ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + EffortIndex {-^ how hard to try -} ->+ ERChebPoly box b ->+ b+chplUpperBoundAffine ix (ERChebPoly coeffs) =+ affiBound coeffs+ where+ affiBound coeffs =+ Map.fold (+) constTerm absCoeffs+ where+ absCoeffs = Map.map abs $ Map.delete chplConstTermKey coeffs+ constTerm = Map.findWithDefault 0 chplConstTermKey coeffs+++{-|+ Find a close upper bound on an affine polynomial over the + unit domain [-1,1]^n. +-}+chplUpperBoundAffineCorners ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box,+ DomainBoxMappable boxb boxbb varid b [(b,b)], Num varid, Enum varid) => + EffortIndex {-^ how hard to try -} ->+ ERChebPoly box b ->+ b+chplUpperBoundAffineCorners ix p@(ERChebPoly coeffs) =+ affiBound (coeffs, vars)+ where+ vars = chplGetVars p+ affiBound (coeffs, vars)+ | null vars =+ Map.findWithDefault 0 chplConstTermKey coeffs+ | otherwise =+ foldl1 max cornerValues+ where+ cornerValues =+ map (\pt -> chplEvalUp pt p) corners+ where+-- corners :: [boxb]+ corners = + map (DBox.fromList . (zip [1..n])) $ prod n+ where+ n = fromInteger $ toInteger $ length vars+ -- n-fold product list of [-1,1]+ prod n + | n == 1 = [[-1],[1]]+ | otherwise =+ (map ((-1):) prodNm1) ++ (map (1:) $ prodNm1)+ where+ prodNm1 = prod (n-1)++{-|+ Find a close upper bound on a quadratic polynomial over the + unit domain [-1,1]^n. +-}+chplUpperBoundQuadr ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box,+ DomainBoxMappable boxra boxras varid (ERInterval b) [ERInterval b],+ DomainBoxMappable boxra boxra varid (ERInterval b) (ERInterval b), + DomainIntBox boxra varid (ERInterval b), Num varid, Enum varid) => + EffortIndex {-^ how hard to try looking for peaks -} ->+ ERChebPoly box b ->+ b+chplUpperBoundQuadr ix p@(ERChebPoly coeffs) =+ quadBound (coeffs, vars)+ where+ vars = chplGetVars p+ quadBound (coeffs, vars)+ | null vars =+ Map.findWithDefault 0 chplConstTermKey coeffs+ | hasInteriorPeak =+ foldl max peakValue edgeBounds+ | otherwise =+ foldl1 max edgeBounds+ where+ edgeBounds =+ map quadBound $ concat $ map removeVar vars+ (hasInteriorPeak, peakValue) =+ case maybePeak of+ Just peak ->+ (noPositiveSquare -- if any term x^2 has a positive coeff, there is no peak + &&+ (and $ map maybeInUnit $ DBox.elems peak)+ ,+ erintv_right $+ chplEvalApprox makeInterval peak p + )+ Nothing -> (False, undefined)+ where+ noPositiveSquare =+ and $ map (<= 0) $ map getQuadCoeff vars+ getQuadCoeff var = + Map.findWithDefault 0 (DBox.singleton var 2) coeffs+ maybeInUnit r =+ case (RA.compareReals r (-1), RA.compareReals (1) r) of+ (Just LT, _) -> False -- ie r < -1+ (_, Just LT) -> False -- ie r > 1+ _ -> True+ maybePeak =+ linearSolver+ (map derivZeroLinearEq vars)+ (DBox.fromList $ map (\v -> (v,(-1) RA.\/ 1)) vars)+ (2^^(-ix))+ where+ derivZeroLinearEq var =+ (linCoeffs, - constCoeff)+ where+ constCoeff =+ makeInterval $+ Map.findWithDefault 0 (DBox.singleton var 1) coeffs+ -- recall T_1(x) = x, T_1'(x) = 1+ linCoeffs =+ DBox.fromList $+ (var, 4 * quadCoeff) -- T_2(x) = 2*x^2 - 1; T_2'(x) = 4*x+ : (map getVarVarCoeff $ var `delete` vars)+ quadCoeff =+ makeInterval $+ Map.findWithDefault 0 (DBox.singleton var 2) coeffs+ getVarVarCoeff var2 =+ (var2,+ makeInterval $+ Map.findWithDefault 0 (DBox.fromList [(var,1), (var2,1)]) coeffs)+ makeInterval b = ERInterval b b+ removeVar var =+ [(substVar True, newVars), + (substVar False, newVars)]+ where+ newVars = var `delete` vars+ substVar isOne =+ chplCoeffs $+ sum $ + map (makeMonomial isOne) $ + Map.toList coeffs+ makeMonomial isOne (term, coeff) =+ ERChebPoly $ Map.fromList $+ case (DBox.toList term) of+ [(v,2)] | v == var ->+ [(chplConstTermKey, coeff)]+ [(v,1)] | v == var ->+ [(chplConstTermKey, + case isOne of True -> coeff; False -> - coeff)]+ [(v1,1), (v2,1)] | v1 == var ->+ [(DBox.fromList [(v2,1)], + case isOne of True -> coeff; False -> - coeff)]+ [(v1,1), (v2,1)] | v2 == var ->+ [(DBox.fromList [(v1,1)], + case isOne of True -> coeff; False -> - coeff)]+ _ ->+ [(term, coeff)]++{-|+ Approximate from below and from above the pointwise maximum of two polynomials+-}+chplMax ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + ERChebPoly box b ->+ ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplMax maxDegree p1 p2 =+ (- (-p1 - differenceDown), p1 + differenceUp)+ where+ (differenceDown, differenceUp) = chplNonneg maxDegree $ p2 - p1++{-|+ Approximate the function max(0,p(x)) by a polynomial from below+ and from above. +-}+chplNonneg ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplNonneg = chplNonnegCubic++{-|+ A version of 'chplNonneg' using a cubic approximation. +-}+chplNonnegCubic ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplNonnegCubic maxDegree p+ | upperB <= 0 = (chplConst 0, chplConst 0)+ | lowerB >= 0 = (p, p)+ | otherwise = -- ie lowerB < 0 < upperB: polynomial may be crossing 0...+ -- work out the cubic polynomial (a3*x^3 + a2*x^2 + a1*x + a0) / b + -- that hits 0 at lowerB with derivative 0 + -- and hits upperB at upperB with derivative 1 + (cubicAppliedOnPDown - valueAt0, cubicAppliedOnPUp + (chplConst correction))+ where + upperB = chplUpperBoundAffine 10 p + lowerB = - (chplUpperBoundAffine 10 (- p))+ cubicAppliedOnPUp = evalCubic multiplyByPUp+ cubicAppliedOnPDown = evalCubic multiplyByPDown+ evalCubic multiplyByP =+ p0 * (chplConst $ recip b)+ where+ p0 = multiplyByP p1 + (chplConst a0) -- ie p*(p*(p * a3 + a2) + a1) + a0+ p1 = multiplyByP p2 + (chplConst a1) -- ie p*(p * a3 + a2) + a1+ p2 = multiplyByP p3 + (chplConst a2) -- ie p * a3 + a2+ p3 = chplConst a3+ multiplyByPUp =+ snd . chplReduceDegree maxDegree . (p *)+ multiplyByPDown =+ fst . chplReduceDegree maxDegree . (p *)+ {-+ The cubic polynomial's coefficients are calculated by solving a system of 4 linear eqs.+ The generic solution is as follows:+ b = (r - l)^3+ a3 = -(r + l)+ a2 = 2*(r^2 + r*l + l^2)+ a1 = -l*(4*r^2 + r*l + l^2)+ a0 = 2*r^2*l^2+ -}+ r = upperB+ l = lowerB+ b = - ((r - l) * ((r - l) * (l - r))) + -- this one has to round downwards because it is a denominator+ a3 = (- r) + (- l) -- remember to round upwards!+ a2 = 2*(r2rll2Up)+ a1 = (- l) * (r2rll2Up + 3*rSqUp) -- since l < 0, the other argument is rounded upwards+ a0 = 2 * rSqUp * lSqUp+ r2rll2Up = rSqUp + r*l + lSqUp + rSqUp = r*r+ lSqUp = l*l+ rSqDown = -((-r)*r)+ lSqDown = -((-l)*l)+ {- + The cubic polynomial may sometimes fail to dominate+ x or sometimes it dips below 0.+ Work out the amount by which it has to be lifted up+ to fix these problems. + -}+ correction+ | 2*rSqDown < l*(r + l) =+ erintv_right $+ (peak0 * (peak0 * (peak0 * (-a3I) - a2I) - a1I) - a0I) / bI+ | 2*lSqDown < r*(r + l) =+ erintv_right $+ ((peakP * (peakP * (peakP * (-a3I) - a2I) - a1I) - a0I) / bI) + peakP+ | otherwise = 0+ where+ -- these have to be computed interval-based:+ [a0I, a1I, a2I, a3I, bI, lI, rI] = + map (\x -> ERInterval x x) [a0,a1,a2,a3,b,l,r]+ peak0 = (lI + 4*rI*rI/(lI+rI)) / 3 + peakP = (rI + 4*lI*lI/(lI+rI)) / 3+ {-+ The same cubic polynomial can be used as a lower bound when+ we subtract its value at 0 rounded upwards.+ -}+ valueAt0 = chplConst $ a0 / b
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs view
@@ -0,0 +1,455 @@+{-# LANGUAGE FlexibleContexts #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary+ Description : (internal) elementary functions applied to polynomials + Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Internal module for "Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom".+ + Implementation of elementary functions applied to polynomials.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Elementary +where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds++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.Base as B+import Data.Number.ER.Real.Approx.Interval+import Data.Number.ER.Real.Arithmetic.Elementary+import qualified Data.Number.ER.Real.DomainBox as DBox+import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainIntBox)+import Data.Number.ER.BasicTypes+import Data.Number.ER.Misc++import qualified Data.Map as Map++{-|+ Approximate the pointwise square root of a polynomial + by another polynomial from below and from above. +-}+chplSqrt ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + EffortIndex {-^ ?? -} -> + ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplSqrt maxDegree ix p =+ error "ERChebPoly: chplSqrt: not implemented yet"++{-|+ Approximate the pointwise exponential of a polynomial + by another polynomial from below and from above. +-}+chplExp ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + EffortIndex {-^ minimum approx Taylor degree -} -> + ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplExp maxDegree ix p =+ (expDownwards, expUpwards + valueAtRDnNeg + (chplConst expRUp))+ where+ expUpwards =+ (chplConst expMUp) * (chplPow maxDegree (expNear0Up pNear0Up) a_int) + expDownwards =+ (chplConst expMDn) * (chplPow maxDegree (expNear0Dn pNear0Dn) a_int) + upperB = chplUpperBoundAffine ix p + lowerB = - (chplUpperBoundAffine ix (- p))+ m = (upperB + lowerB) / 2+ r = (upperB - lowerB) / 2 + expMUp = erintv_right expM + expMDn = erintv_left expM+ expM = erExp_R ix (ERInterval m m)+ pNear0Up = (p - (chplConst m)) * (chplConst $ recip a_base)+ pNear0Dn = - (((-p) + (chplConst m)) * (chplConst $ recip a_base))+ a_base = fromInteger a_int+ a_int = max 1 $ floor r -- could this be too high?+ expNear0Up p0 =+ expAux p0 1 (B.setGranularity coeffGr 1)+ expNear0Dn p0 =+ negate $ expAux p0 1 (B.setGranularity coeffGr (-1))+ expAux p0 nextDegree thisCoeff+ | nextDegree > taylorDegree =+ chplConst thisCoeff+ | otherwise =+ snd $ chplReduceDegree maxDegree $+ (chplConst thisCoeff) + p0 * (expAux p0 (nextDegree + 1) nextCoeff)+ where+ nextCoeff = + thisCoeff / (fromInteger nextDegree)+ taylorDegree = 1 + 2 * (ix `div` 6)+ coeffGr = effIx2gran $ 10 + 3 * taylorDegree+ expRUp = erintv_right expR+ expR = erExp_R ix (ERInterval r r)+ valueAtRDnNeg = + expAux (chplConst r) 1 (B.setGranularity coeffGr (-1))++ +{-|+ Approximate the pointwise integer power of a polynomial by another polynomial from above. +-}+chplPow ::+ (B.ERRealBase b, Integral i, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + ERChebPoly box b ->+ i ->+ ERChebPoly box b+chplPow maxDegree p n+ | n == 0 =+ chplConst 1+ | n == 1 =+ p + | even n =+ snd $ chplReduceDegree maxDegree $ powHalfN * powHalfN+ | odd n =+ snd $ chplReduceDegree maxDegree $ + p * + (snd $ chplReduceDegree maxDegree $+ powHalfN * powHalfN)+ where+ powHalfN =+ chplPow maxDegree p halfN+ halfN = n `div` 2+ +{-|+ Approximate the pointwise natural logarithm of a polynomial + by another polynomial from below and from above. +-}+chplLog ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + EffortIndex {-^ ?? -} -> + ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplLog maxDegree ix p =+ error "ERChebPoly: chplLog: not implemented yet"++{-|+ Approximate the pointwise sine of a polynomial + by another polynomial from below and from above. +-}+chplSineCosine ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) =>+ Bool {-^ True iff sine, False iff cosine -} -> + Int {-^ maximum polynomial degree -} -> + EffortIndex {-^ minimum approx Taylor degree -} -> + ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplSineCosine isSine maxDegree ix p+ -- p - 2k*pi range within [-pi/2, pi/2]: + | ranfNear0 `RA.refines` plusMinusPiHalf =+-- unsafePrint+-- (+-- "ERChebPoly: chplSineCosine: [-pi/2, pi/2]: "+-- ++ "\n p = " ++ show p+-- ++ "\n ranf = " ++ show ranf+-- ++ "\n k = " ++ show k+-- ++ "\n ranfNear0 = " ++ show ranfNear0+-- ) $+ case isSine of+ True -> sineShifted (- k2pi)+ False -> cosineShifted (- k2pi)+ -- p - 2k*pi range within [0, pi]: + | (ranfNear0 - piHalf) `RA.refines` plusMinusPiHalf =+-- unsafePrint+-- (+-- "ERChebPoly: chplSineCosine: [0, pi]: "+-- ++ "\n p = " ++ show p+-- ++ "\n ranf = " ++ show ranf+-- ++ "\n k = " ++ show k+-- ++ "\n ranfNear0 = " ++ show ranfNear0+-- ) $+ case isSine of+ -- use sin(x) = cos(x - pi/2) and cos(x) = - sin(x - pi/2):+ True -> cosineShifted (- k2pi - piHalf)+ False -> sineShiftedNegated (- k2pi - piHalf)+ -- p - 2k*pi range within [-pi, 0]: + | (ranfNear0 + piHalf) `RA.refines` plusMinusPiHalf =+ case isSine of+ -- use sin(x) = - cos(x + pi/2) and cos(x) = sin(x + pi/2):+ True -> cosineShiftedNegated (-k2pi + piHalf)+ False -> sineShifted (-k2pi + piHalf)+ -- p - 2k*pi range within [pi/2, 3pi/2]: + | (ranfNear0 - pi) `RA.refines` plusMinusPiHalf =+ -- use sin(x) = - sin(x - pi) and cos(x) = - cos(x - pi)+ case isSine of+ True -> sineShiftedNegated (- k2pi - pi)+ False -> cosineShiftedNegated (- k2pi - pi)+ | otherwise = (chplConst (-1), chplConst 1)+-- (expDownwards, expUpwards + valueAtRDnNeg + (chplConst expRUp))+ where+ ranfNear0 = ranf - k2pi + k2pi = k * 2 * pi+ plusMinusPiHalf = (-piHalfLO) RA.\/ piHalfLO+ pi = RAEL.pi ix + piHalf = pi / 2+ (piHalfLO, piHalfHI) = RA.bounds piHalf+ ranf = + ERInterval + (negate $ chplUpperBoundAffine 10 (-p)) + (chplUpperBoundAffine 10 p)+ k = + fromInteger $ floor $ + case (pi + ranf) / (2 * pi) of ERInterval lo hi -> lo+ + sineShiftedNegated shift =+ boundsNegate $ sineShifted shift+ + cosineShiftedNegated shift =+ boundsNegate $ cosineShifted shift++ boundsNegate (pLO, pHI) = (- pHI, - pLO)+ + sineShifted shift =+ boundsAddErr shiftWidthB $ sineTaylor (p + shiftPoly) (ranf + shift)+ where+ shiftPoly = chplConst shiftLOB+ ERInterval shiftLOB shiftHIB = shift+ shiftWidthB = shiftHIB - shiftLOB+ + cosineShifted shift =+ boundsAddErr shiftWidthB $ cosineTaylor (p + shiftPoly) (ranf + shift)+ where+ shiftPoly = chplConst shiftLOB+ ERInterval shiftLOB shiftHIB = shift+ shiftWidthB = shiftHIB - shiftLOB+ + boundsAddErr errB (pLO, pHI) =+ (pLO `plusDown` (- errPoly), pHI + errPoly)+ where+ errPoly = chplConst errB+ + sineTaylor x xran =+ (sineDown, sineUp)+ where+ sineUp =+ chplReduceDegreeUp maxDegree $ + x * sineUpTaylor + (chplConst sineUpErrorBound)+ (sineUpTaylor, sineUpErrorTermDegree, sineUpErrorTermCoeff) =+ taylorAux x 1 (B.setGranularity coeffGr 1)+ sineUpErrorBound =+ case sineUpErrorBoundRA of ERInterval lo hi -> hi+ where+ sineUpErrorBoundRA = + (xranLargerEndpoint ^ (1 + sineUpErrorTermDegree)) * sineUpErrorTermCoeffRA+ sineUpErrorTermCoeffRA =+ abs $+ ERInterval sineUpErrorTermCoeff sineUpErrorTermCoeff+ sineDown = + negate $ chplReduceDegreeUp maxDegree $ + x * sineDownTaylorNeg + (chplConst $ sineDownErrorBound)+ (sineDownTaylorNeg, sineDownErrorTermDegree, sineDownErrorTermCoeff) =+ taylorAux x 1 (B.setGranularity coeffGr (-1))+ sineDownErrorBound =+ case sineDownErrorBoundRA of ERInterval lo hi -> hi+ where+ sineDownErrorBoundRA =+ (xranLargerEndpoint ^ (1 + sineDownErrorTermDegree)) * sineDownErrorTermCoeffRA+ sineDownErrorTermCoeffRA =+ abs $+ ERInterval sineDownErrorTermCoeff sineDownErrorTermCoeff+ xranLargerEndpoint = + max (abs xranLO) (abs xranHI)+ (xranLO, xranHI) = RA.bounds xran+ + cosineTaylor x xran =+-- unsafePrint+-- (+-- "ERChebPoly.Elementary: chplSineCosine: cosineTaylor: "+-- ++ "\n xran = " ++ show xran+-- ++ "\n cosineUpErrorBound = " ++ show cosineUpErrorBound+-- ++ "\n cosineUpErrorTermDegree = " ++ show cosineUpErrorTermDegree+-- ++ "\n cosineUpErrorTermCoeff = " ++ show cosineUpErrorTermCoeff+-- ++ "\n xranLargerEndpoint = " ++ show xranLargerEndpoint+-- )+ (cosineDown, cosineUp)+ where+ cosineUp =+ chplReduceDegreeUp maxDegree $ + cosineUpTaylor + (chplConst cosineUpErrorBound)+ (cosineUpTaylor, cosineUpErrorTermDegree, cosineUpErrorTermCoeff) =+ taylorAux x 0 (B.setGranularity coeffGr 1)+ cosineUpErrorBound+ | odd (cosineUpErrorTermDegree `div` 2) = 0+ | otherwise =+ case cosineUpErrorBoundRA of ERInterval lo hi -> hi+ where+ cosineUpErrorBoundRA = + (xranLargerEndpoint ^ (cosineUpErrorTermDegree)) * cosineUpErrorTermCoeffRA+ cosineUpErrorTermCoeffRA =+ abs $+ ERInterval cosineUpErrorTermCoeff cosineUpErrorTermCoeff+ cosineDown = + negate $ chplReduceDegreeUp maxDegree $ + cosineDownTaylorNeg + (chplConst $ cosineDownErrorBound)+ (cosineDownTaylorNeg, cosineDownErrorTermDegree, cosineDownErrorTermCoeff) =+ taylorAux x 0 (B.setGranularity coeffGr (-1))+ cosineDownErrorBound + | even (cosineDownErrorTermDegree `div` 2) = 0+ | otherwise =+ case cosineDownErrorBoundRA of ERInterval lo hi -> hi+ where+ cosineDownErrorBoundRA =+ (xranLargerEndpoint ^ (cosineDownErrorTermDegree)) * cosineDownErrorTermCoeffRA+ cosineDownErrorTermCoeffRA =+ abs $+ ERInterval cosineDownErrorTermCoeff cosineDownErrorTermCoeff+ xranLargerEndpoint = + max (abs xranLO) (abs xranHI)+ (xranLO, xranHI) = RA.bounds xran+ + taylorAux p0 thisDegree thisCoeff+ | nextDegree > taylorDegree =+-- unsafePrint+-- (+-- "ERChebPoly: chplSine: taylorAux: "+-- ++ "\n thisCoeff = " ++ show thisCoeff+-- ++ "\n nextDegree = " ++ show nextDegree+-- )+ (chplConst thisCoeff, nextDegree, nextCoeff)+ | otherwise =+-- unsafePrint+-- (+-- "ERChebPoly: chplSine: taylorAux: "+-- ++ "\n thisCoeff = " ++ show thisCoeff+-- ++ "\n nextDegree = " ++ show nextDegree+-- ++ "\n errorTermCoeff = " ++ show errorTermCoeff+-- ++ "\n errorTermDegree = " ++ show errorTermDegree+-- )+ (chplReduceDegreeUp maxDegree $+ (chplConst thisCoeff) + p0 * p0 * rest,+ errorTermDegree, errorTermCoeff) + where+ (rest, errorTermDegree, errorTermCoeff) =+ taylorAux p0 nextDegree nextCoeff+ nextDegree = thisDegree + 2+ nextCoeff = + thisCoeff / (fromInteger $ negate $ nextDegree * (nextDegree - 1))+ taylorDegree = ix `div` 3+ coeffGr = effIx2gran $ ix++{-|+ Approximate the pointwise cosine of a polynomial + by another polynomial from below and from above+ using the tau method + as described in [Mason & Handscomb 2003, p 62]. +-}+chplRecip ::+ (B.ERRealBase b, RealFrac b, DomainBox box varid Int, Ord box) => + Int {-^ maximum polynomial degree -} -> + EffortIndex {-^ minimum approx degree -} -> + ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplRecip maxDegree ix p+ | upperB < 0 = -- range negative+ (\(lo, hi) -> (-hi, -lo)) $ chplRecip maxDegree ix (negate p)+ | lowerB > 0 = -- range positive+-- unsafePrint+-- (+-- "ERChebPoly: chplRecip: "+-- ++ "\n k = " ++ show k+-- ++ "\n lowerB = " ++ show lowerB+-- ++ "\n tau = " ++ (show $ recip tauInv)+-- )+ (resDn, resUp)+ | otherwise = -- cannot establish 0 freedom+ error $+ "ERChebPoly: chplRecip: "+ ++ "cannot deal with estimated range " ++ show ranp+ ++ "of polynomial: \n" ++ show p + where+ ranp = ERInterval lowerB upperB+ lowerB = - (chplUpperBoundAffine ix (- p))+ upperB = chplUpperBoundAffine ix p+ + tauDegree = effIx2int (ix `div` 3)+ coeffGr = effIx2gran $ ix+ + -- translate p to have range above 1:+ k = intLogUp 2 $ ceiling (1/lowerB) -- 2^k * lowerB >= 1+ upperBtr = upperB * 2^k -- upper bound of translated poly+ (pAbove1Dn, pAbove1Up) = -- p multiplied by 2^k; range in [1,upperBtr] + chplScale (2^k) p+ + -- translate T_1 to domain [0, upperBtr] and apply it to (pAbove1 - 1):+ -- T'_1 = nu * (p - 1) - 1+ trT1Dn = + (chplScaleDown nuLOB (pAbove1Dn - 1)) - 1+ trT1Up =+ (chplScaleUp nuHIB (pAbove1Up - 1)) - 1+ nu = recip nuInv -- auxiliary constant+ ERInterval nuLOB nuHIB = nu+ nuInv = (RA.setMinGranularity coeffGr (ERInterval upperBtr upperBtr)) / 2+ nuPlus1 = nu + 1+ nuInvPlus1 = nuInv + 1+ nuInvDiv2 = nuInv / 2+ + -- define such translated T_i's for all i >= 0:+ trTis =+ map (mapPair (chplReduceDegreeDown maxDegree, chplReduceDegreeUp maxDegree)) $ + chebyEvalTsRoundDownUp trT1Dn + + -- construct the result from interval coefficients:+ resDn =+ chplScaleDown (2^k) $+ (-tauAbsUpPoly) `plusDown` + (chplScaleUp tauAbsDnB $+ sumDown $+ (- errPoly) : (zipWith scaleDn cis trTis))+ resUp =+ chplScaleUp (2^k) $+ (tauAbsUpPoly) `plusUp` + (chplScaleUp tauAbsUpB $+ sumUp $+ (errPoly) : (zipWith scaleUp cis trTis))+ + scaleDn c (trTDn, trTUp) + | r >= 0 = chplScaleDown r trTDn+ | otherwise = chplScaleDown r trTUp+ where+ r = c * tauSign+ scaleUp c (trTDn, trTUp) + | r >= 0 = chplScaleUp r trTUp+ | otherwise = chplScaleUp r trTDn+ where+ r = c * tauSign+ + tauAbsUpPoly = chplConst $ tauAbsUpB+ tauSign = + case RA.compareReals tauInv 0 of+ Just GT -> 1+ Just LT -> -1+ ERInterval tauAbsDnB tauAbsUpB = abs $ recip tauInv+ cis =+ map (\(ERInterval lo hi) -> hi) c0n + errPoly = chplConst err+ err =+ foldl1 plusUp $+ map (\(ERInterval lo hi) -> hi - lo) c0n+ + -- work out the coefficients in interval arithmetic using the tau method:+ c0n = c0 : c1n+ tauInv = c0 * nuInvPlus1 + c1 * nuInvDiv2+ c0 = - c1 * nuPlus1 - c2/2+ (c1 : c2 : _) = c1n+ c1n = reverse $ take n $ csRev+ n = tauDegree+ csRev =+ cn : cnM1 : (csAux cn cnM1)+ where+ cn = 1+ cnM1 = - 2 * nuPlus1+ csAux cn cnM1 =+ cnM2 : (csAux cnM1 cnM2)+ where+ cnM2 = - cn - 2 * nuPlus1 * cnM1
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs view
@@ -0,0 +1,190 @@+{-# LANGUAGE FlexibleContexts #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval+ Description : (internal) evaluation of polynomials+ Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++ Internal module for "Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom".+ + Implementation of various evaluation functions related to polynomials.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval +where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field++import qualified Data.Number.ER.Real.Approx as RA+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 qualified Data.Map as Map++{-|+ Evaluate a polynomial at a point, consistently rounding upwards and downwards. +-}+chplEval ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box, + DomainBoxMappable boxb boxbb varid b [(b,b)]) => + boxb -> + ERChebPoly box b ->+ (b, b)+chplEval vals (ERChebPoly coeffs) =+ (foldl plusDown 0 termValsLo, foldl plusUp 0 termValsHi)+ where+ (termValsLo, termValsHi) =+ unzip $ map evalTerm $ Map.toList coeffs+ evalTerm (term, c) =+ (foldl timesDown c valsLo, foldl timesUp c valsHi)+ where+ (valsLo, valsHi) = + unzip $ map evalVar $ DBox.toList term+ evalVar (varID, degree) =+ (DBox.lookup "ERChebPoly.Eval: chplEval" varID valsDegrees) !! degree+ valsDegrees =+ DBox.map chebyEvalTsRoundDownUp vals++chplEvalDown, chplEvalUp ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box, + DomainBoxMappable boxb boxbb varid b [(b,b)]) => + boxb -> + ERChebPoly box b ->+ b+chplEvalUp pt = snd . chplEval pt+chplEvalDown pt = fst . chplEval pt++chebyEvalTsRoundDownUp ::+ (Num v) =>+ v -> [(v,v)]+chebyEvalTsRoundDownUp val =+ chebyIterate (1,1) (val, val)+ where+ chebyIterate tNm2@(tNm2Down, tNm2Up) tNm1@(tNm1Down, tNm1Up) =+ tNm2 : (chebyIterate tNm1 (tNDown, tNUp))+ where+ tNUp = 2 * val * tNm1Up - tNm2Down + tNDown = ((2 * val) `timesDown` tNm1Down) - tNm2Up ++chebyEvalTsExact ::+ (Num v) =>+ v -> [v] +chebyEvalTsExact val =+ chebyIterate 1 val+ where+ chebyIterate tNm2 tNm1 =+ tNm2 : (chebyIterate tNm1 tN)+ where+ tN = 2 * val * tNm1 - tNm2 ++{-|+ Evaluate a polynomial at a real number approximation +-}+chplEvalApprox ::+ (B.ERRealBase b, RA.ERApprox ra, + DomainBox box varid Int, Ord box,+ DomainBoxMappable boxra boxras varid ra [ra], + DomainIntBox boxra varid ra) =>+ (b -> ra) -> + boxra -> + ERChebPoly box b ->+ ra+chplEvalApprox b2ra vals (ERChebPoly coeffs) =+ sum $ map evalTerm $ Map.toList coeffs+ where+ evalTerm (term, c) =+ (b2ra c) * (product $ map evalVar $ DBox.toList term)+ evalVar (varID, degree) =+ (DBox.lookup "ERChebPoly.Eval: chplEvalApprox: " varID valsDegrees) !! degree+ valsDegrees =+ DBox.map chebyEvalTsExact vals+ +{-|+ Substitute several variables in a polynomial with real number approximations,+ rounding downwards and upwards.+-}+chplPartialEvalApprox ::+ (B.ERRealBase b, RA.ERApprox ra, + DomainBox box varid Int, Ord box,+ DomainBoxMappable boxra boxras varid ra [ra], + DomainIntBox boxra varid ra) =>+ (ra -> (b,b)) ->+ boxra ->+ ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplPartialEvalApprox ra2endpts substitutions (ERChebPoly coeffs) =+ (ERChebPoly $ Map.insertWith plusDown chplConstTermKey (- corr) coeffsSubstDown, + ERChebPoly $ Map.insertWith plusUp chplConstTermKey corr coeffsSubstUp)+ where+ (coeffsSubstDown, coeffsSubstUp, corr) =+ Map.foldWithKey processTerm (Map.empty, Map.empty, 0) coeffs+ processTerm termKey coeff (coeffsSubstDownPrev, coeffsSubstUpPrev, corrPrev) =+ (Map.insertWith plusDown newTermKey newCoeffDown coeffsSubstDownPrev,+ Map.insertWith plusUp newTermKey newCoeffUp coeffsSubstUpPrev,+ corrPrev + corrVars)+ where+ corrVars = (substValHi - substValLo) * coeff+ (newCoeffDown, newCoeffUp) + | coeff > 0 = (coeff `timesDown` substValLo, coeff `timesUp` substValHi)+ | coeff < 0 = (coeff `timesDown` substValHi, coeff `timesUp` substValLo)+ | otherwise = (0,0)+ (substValLo, substValHi) = ra2endpts substVal+ (substVal, newTermKey) =+ DBox.foldWithKey processVar (1, DBox.noinfo) termKey+ processVar varID degree (substValPrev, newTermKeyPrev) =+ case DBox.member varID substitutions of+ True -> + (substValPrev * (evalVar varID degree), newTermKeyPrev)+ False ->+ (substValPrev, DBox.insert varID degree newTermKeyPrev)+ evalVar varID degree =+ (DBox.lookup "ERChebPoly.Eval: chplPartialEvalApprox: " varID valsDegrees) !! degree+ valsDegrees =+ DBox.map chebyEvalTsExact substitutions+ ++{-|+ Compose two polynomials, rounding upwards+ provided the second polynomial maps [-1,1] into [-1,1].+-}+chplCompose ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ Int ->+ ERChebPoly box b ->+ Map.Map varid (ERChebPoly box b) + {-^ variable to substitute, polynomial to substitute -} ->+ (ERChebPoly box b, ERChebPoly box b)+chplCompose maxDegree p@(ERChebPoly coeffs) substitutions =+ (foldl plusDown 0 termValsLo, foldl plusUp 0 termValsHi)+ where+ (termValsLo, termValsHi) =+ unzip $ map evalTerm $ Map.toList coeffs+ evalTerm (term, c) =+ (foldl timesDown cPoly valsLo, foldl timesUp cPoly valsHi)+ where+ cPoly = chplConst c+ (valsLo, valsHi) = + unzip $ map evalVar $ DBox.toList term+ evalVar (varID, degree) =+ case Map.lookup varID substDegrees of+ Nothing ->+ (varPoly, varPoly)+ Just pvDegrees ->+ pvDegrees !! degree+ where+ varPoly = + ERChebPoly $ Map.singleton (DBox.singleton varID degree) 1+ substDegrees =+ Map.map mkPVDegrees substitutions+ mkPVDegrees pv =+ map + (mapPair + (chplReduceDegreeDown maxDegree, + chplReduceDegreeUp maxDegree)) $ + chebyEvalTsRoundDownUp pv
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs view
@@ -0,0 +1,226 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field+ Description : (internal) field operations applied to polynomials + Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Internal module for "Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom".+ + Implementation of field arithmetic over polynomials + with rounding consistent over the whole domain.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field ++where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic++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, DomainIntBox)+import Data.Number.ER.Misc++import qualified Data.Map as Map++chplAffine ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ b -> + Map.Map varid b ->+ ERChebPoly box b+chplAffine at0 varCoeffs =+ ERChebPoly $ + Map.insert chplConstTermKey at0 $+ Map.mapKeys (\ i -> DBox.singleton i 1) varCoeffs+ +{-|+ Convert a polynomial to a lower-order one that is dominated by (resp. dominates)+ it closely on the domain [-1,1].+-}+chplReduceDegree ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + Int {-^ new maximal order -} ->+ ERChebPoly box b -> + (ERChebPoly box b, ERChebPoly box b) {-^ lower and upper bounds with limited degree -}+chplReduceDegree maxOrder (ERChebPoly coeffs) =+ (ERChebPoly newCoeffsDown, ERChebPoly newCoeffsUp)+-- errorModule "chplSetMaxOrder: not implemented yet"+ where+ newCoeffsUp =+ Map.insertWith plusUp chplConstTermKey highOrderCompensation coeffsLowOrder+ newCoeffsDown =+ Map.insertWith plusDown chplConstTermKey (-highOrderCompensation) coeffsLowOrder+ highOrderCompensation =+ Map.fold (\ new prev -> prev + (abs new)) 0 coeffsHighOrder+ (coeffsHighOrder, coeffsLowOrder) = + Map.partitionWithKey (\ k v -> chplTermOrder k > maxOrder) coeffs++chplReduceDegreeDown m = fst . chplReduceDegree m+chplReduceDegreeUp m = snd . chplReduceDegree m++instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Num (ERChebPoly box b)+ where+ fromInteger n =+ ERChebPoly $ Map.singleton chplConstTermKey (fromInteger n)+ abs (ERChebPoly coeffs) =+ errorModule "abs of a polynomial not implemented, use UFB.max instead"+ signum (ERChebPoly coeffs) =+ errorModule "signum of a polynomial not implemented, use RA.leqReals instead"+ --------- negation ----------+ negate (ERChebPoly coeffs) =+ ERChebPoly $ Map.map negate coeffs+ --------- addition ----------+ (ERChebPoly coeffs1) + (ERChebPoly coeffs2) =+ ERChebPoly sumCoeffs+ where+ sumCoeffs =+ Map.insertWith (+) chplConstTermKey maxError coeffsDown+ -- point-wise sum of polynomials with coeff rounding errors+ -- compensated for by enlarging the constant term+ coeffsUp =+ (Map.unionWith (+) coeffs1 coeffs2)+ -- point-wise sum of polynomials with coeffs rounded upwards+ coeffsDown =+ (Map.unionWith (\c1 c2 -> - ((- c1) + (- c2))) coeffs1 coeffs2)+ -- point-wise sum of polynomials with coeffs rounded upwards+ maxError =+ Map.fold (+) 0 $ + Map.intersectionWith (-) coeffsUp coeffsDown+ -- addition must round upwards on interval [-1,1]+ -- non-constant terms are multiplied by quantities in [-1,1] + -- and thus can make the result drop below the exact result+ -- -> to compensate add the rounding difference to the constant term + --------- multiplication ----------+ (ERChebPoly coeffs1) * (ERChebPoly coeffs2) =+ ERChebPoly prodCoeffs+ where + prodCoeffs =+ Map.insertWith (+) chplConstTermKey roundOffCompensation $ + Map.map negate directProdCoeffsDown+ roundOffCompensation =+ Map.fold (+) 0 $+ Map.unionWith (+) directProdCoeffsDown directProdCoeffsUp+ (directProdCoeffsUp, directProdCoeffsDown) =+ foldl addCombiCoeff (Map.empty, Map.empty) combinedCoeffs+ where+ addCombiCoeff+ (prevCoeffsUp, prevCoeffsDown) + (coeffUp, coeffDown, (powersList, coeffCount)) =+ foldl addOnce (prevCoeffsUp, prevCoeffsDown) powersList+ where+ addOnce (prevCoeffsUp, prevCoeffsDown) powers =+ (Map.insertWith (+) powers coeffUpFrac prevCoeffsUp, + Map.insertWith (+) powers coeffDownFrac prevCoeffsDown)+ coeffUpFrac = coeffUp / coeffCountB+ coeffDownFrac = coeffDown / coeffCountB+ coeffCountB = fromInteger coeffCount+ combinedCoeffs =+ [ -- (list of triples)+ (+ (c1 * c2) -- upwards rounded product+ ,+ ((- c1) * c2) -- downwards rounded negated product+ ,+ combinePowers powers1 powers2+ )+ |+ (powers1, c1) <- coeffs1List,+ (powers2, c2) <- coeffs2List+ ]+ combinePowers powers1 powers2 =+ (combinedPowers, 2 ^ (length sumsDiffs)) + where+ combinedPowers =+ map (DBox.fromAscList . (filter $ \ (k,v) -> v > 0)) $+ allPairsCombinations $ + sumsDiffs+ sumsDiffs = + -- associative list with the sum and difference of powers for each variable+ zipWith (\(k,s) (_,d) -> (k,(s,d)))+ (DBox.toAscList $ DBox.unionWith (\a b -> (a + b)) powers1 powers2)+ (DBox.toAscList $ DBox.unionWith (\a b -> abs (a - b)) powers1 powers2)+ coeffs1List =+ Map.toList coeffs1+ coeffs2List =+ Map.toList coeffs2+++-- | multiply a polynomial by a scalar rounding downwards and upwards +chplScale ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ b -> + (ERChebPoly box b) -> + (ERChebPoly box b, ERChebPoly box b)+chplScale ratio (ERChebPoly coeffs) =+ (ERChebPoly coeffsDown, ERChebPoly coeffsUp)+ where+ coeffsDown = + Map.insertWith plusDown chplConstTermKey (- errBound) coeffsScaled + coeffsUp = + Map.insertWith plusUp chplConstTermKey errBound coeffsScaled+ (errBound, coeffsScaled) =+ Map.mapAccum processTerm 0 coeffs+ processTerm errBoundPrev coeff =+ (errBoundPrev + errBoundHere, coeffScaledUp)+ where+ errBoundHere = coeffScaledUp - coeffScaledDown+ coeffScaledDown = ratio `timesDown` coeff+ coeffScaledUp = ratio `timesUp` coeff ++chplScaleDown r = fst . chplScale r+chplScaleUp r = snd . chplScale r++-- | multiply a polynomial by a scalar interval rounding downwards and upwards +chplScaleApprox :: + (B.ERRealBase b, DomainBox box varid Int, Ord box) =>+ (b, b) -> + (ERChebPoly box b) -> + (ERChebPoly box b, ERChebPoly box b)+chplScaleApprox (ratioDown, ratioUp) (ERChebPoly coeffs) =+ (ERChebPoly coeffsDown, ERChebPoly coeffsUp)+ where+ coeffsDown =+ Map.insertWith plusDown chplConstTermKey (- errBound) coeffsScaled + coeffsUp = + Map.insertWith plusUp chplConstTermKey errBound coeffsScaled+ (errBound, coeffsScaled) =+ Map.mapAccum processTerm 0 coeffs+ processTerm errBoundPrev coeff =+ (errBoundPrev + errBoundHere, coeffScaledUp)+ where+ errBoundHere = coeffScaledUp - coeffScaledDown+ (coeffScaledDown, coeffScaledUp)+ | coeff >= 0 = + (ratioDown `timesDown` coeff, ratioUp `timesUp` coeff)+ | coeff < 0 = + (ratioUp `timesDown` coeff, ratioDown `timesUp` coeff)+++instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Fractional (ERChebPoly box b)+ where+ fromRational r =+ ERChebPoly $ Map.singleton chplConstTermKey (fromRational r)+ --------- division ----------+ _ / _ =+ errorModule "for division use chplRecip from module Elementary" + +instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Ord (ERChebPoly box b)+ where+ compare _ _ =+ errorModule "cannot compare polynomials, consider using RA.leqReals instead"+ +--instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => Real (ERChebPoly box b)+-- where+-- toRational _ =+-- errorModule "toRational: cannot convert polynomial to rational" +-- +--instance (B.ERRealBase b, DomainBox box varid Int, Ord box) => RealFrac (ERChebPoly box b)+-- where+-- properFraction _ =+-- errorModule "properFraction: rounding of polynomials not implemented" +
+ src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs view
@@ -0,0 +1,169 @@+{-# LANGUAGE FlexibleContexts #-}+{-|+ Module : Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Integration+ Description : (internal) integration of polynomials etc + Copyright : (c) 2007-2008 Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable+ + Internal module for "Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom".+ + Implementation of safely rounded integration of polynomials+ and other related functions.+-}+module Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Integration +where++import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Basic+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Eval+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Field+import Data.Number.ER.RnToRm.UnitDom.ChebyshevBase.Polynom.Bounds++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 qualified Data.Map as Map++{-|+ Approximate from below and from above the integral of a polynomial.+ + Based on the following formulas for Chebyshev polynomials:+ +> \int T_n(x)dx = +> T_{n+1}(x)/2(n+1) - T_{n-1}(x)/2(n-1)++> \int T_1(x)dx = +> T_2(x)/4 + 1/4++> \int T_0(x)dx = +> T_1(x)+ +-}+chplIntegrate ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + varid {-^ variable to integrate by -} -> + ERChebPoly box b ->+ (ERChebPoly box b, ERChebPoly box b)+chplIntegrate x (ERChebPoly coeffs) =+-- unsafePrint+-- (+-- "ERChebPoly: integrate:"+-- ++ "\n pNp1Down = " ++ chplShow True pNp1Down +-- ++ "\n pNm1Up = " ++ chplShow True pNm1Up +-- )+ (chplNormaliseDown $ pNp1Down - pNm1Up, + chplNormaliseUp $ pNp1Up - pNm1Down)+ where+ pNp1Up =+ ERChebPoly $ + Map.insertWith plusUp chplConstTermKey errBoundNp1 $ + Map.fromList coeffsNp1+ pNp1Down =+ ERChebPoly $ + Map.insertWith plusDown chplConstTermKey (- errBoundNp1) $ + Map.fromList coeffsNp1+ pNm1Up =+ ERChebPoly $ + Map.insertWith plusUp chplConstTermKey errBoundNm1 $ + Map.fromList coeffsNm1+ pNm1Down =+ ERChebPoly $ + Map.insertWith plusDown chplConstTermKey (- errBoundNm1) $ + Map.fromList coeffsNm1+ (coeffsNp1, errBoundNp1) =+ foldl cfNp1 ([],0) coeffsList+ (coeffsNm1, errBoundNm1) =+ foldl cfNm1 ([],0) coeffsList+ coeffsList = Map.toList coeffs+ cfNp1 (prevTerms, prevErr) (termKey, coeff)+ | n == 0 =+ ((termKeyNp1, coeff):prevTerms, prevErr)+ | n == 1 =+ ((termKeyNm1, coeff0Up):(termKeyNp1, coeffNp1Up):prevTerms, prevErr + coeff0Err + coeffNp1Err)+ | otherwise =+ ((termKeyNp1, coeffNp1Up):prevTerms, prevErr + coeffNp1Err)+ where+ termKeyNp1 = DBox.insert x (n + 1) termKey+ termKeyNm1 = DBox.insert x (n - 1) termKey + n = DBox.findWithDefault 0 x termKey+ coeffNp1Err = coeffNp1Up - coeffNp1Down + coeffNp1Up = coeff / (2*nB + 2)+ coeffNp1Down = -((-coeff) / (2*nB + 2))+ nB = fromInteger $ toInteger n+ coeff0Up = coeff / 4+ coeff0Down = - ((- coeff) / 4)+ coeff0Err = coeff0Up - coeff0Down + cfNm1 (prevTerms, prevErr) (termKey, coeff)+ | n == 0 || n == 1 =+ ((chplConstTermKey, 0):prevTerms, prevErr)+ | otherwise =+ ((termKeyNm1, coeffNm1Up):prevTerms, prevErr + coeffNm1Err)+ where+ termKeyNm1 = DBox.insert x (n - 1) termKey + n = DBox.findWithDefault 0 x termKey+ coeffNm1Up = coeff / (2*nB - 2)+ coeffNm1Down = -((-coeff) / (2*nB - 2))+ nB = fromInteger $ toInteger n+ coeffNm1Err = coeffNm1Up - coeffNm1Down ++{-|+ measure the volume between a polynomial and the zero axis on [-1,1]^n+-}+chplVolumeAboveZero ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box, + DomainBoxMappable boxb boxbb varid b [(b,b)]) =>+ [varid] ->+ ERChebPoly box b ->+ (b,b)+chplVolumeAboveZero vars p@(ERChebPoly coeffs) =+-- unsafePrint ("chplVolumeAboveZero: returning:" ++ show result) $+-- unsafePrint ("chplVolumeAboveZero: vars = " ++ show vars) $+ result+ where+ result = + (- (integUpAtOddCorners - integDownAtEvenCorners), integUpAtEvenCorners - integDownAtOddCorners)+ integUpAtEvenCorners = sumUp $ map (\pt -> chplEvalUp pt integUp) evenCorners+ integUpAtOddCorners = sumUp $ map (\pt -> chplEvalUp pt integUp) oddCorners + integDownAtEvenCorners = sumDown $ map (\pt -> chplEvalDown pt integDown) evenCorners + integDownAtOddCorners = sumDown $ map (\pt -> chplEvalDown pt integDown) oddCorners+ evenCorners = map (DBox.fromList) evenCornersL+ oddCorners = map (DBox.fromList) oddCornersL+ (evenCornersL, oddCornersL) =+ allPairsCombinationsEvenOdd $ zip vars $ repeat (1,-1)+ integUp = integrateByAllVars snd p vars+ integDown = integrateByAllVars fst p vars+ integrateByAllVars pick p [] = p+ integrateByAllVars pick p (x : xs) =+ integrateByAllVars pick ip xs+ where+ ip = pick $ chplIntegrate x p+-- vars = chplGetVars p+ + +--{-|+-- Calculate approximations to the Chebyshev nodes.+---}+--chebNodes ::+-- (B.ERRealBase b) =>+-- Granularity ->+-- [[b]] -- ^ ith element is the ordered list of ith order Chebyshev nodes +--chebNodes gran =+-- error "ERChebPoly: chebNodes: not implemented yet"+ + +{-|+ Differentiate a polynomial using one of its variables. +-}+chplDifferentiate ::+ (B.ERRealBase b, DomainBox box varid Int, Ord box) => + ERChebPoly box b ->+ varid {-^ variable to differentiate over -} ->+ ERChebPoly box b+chplDifferentiate (ERChebPoly coeffs) varName =+ errorModule "chplDifferentiate: not implemented yet"+
+ tests/Demo.hs view
@@ -0,0 +1,114 @@+{-| + Module : Main+ Description : simple examples of using AERN-RnToRm+ Copyright : (c) Michal Konecny+ License : BSD3++ Maintainer : mik@konecny.aow.cz+ Stability : experimental+ Portability : portable++ Simple examples of using AERN-RnToRm.+-}+module Main where++import qualified Data.Number.ER.RnToRm as AERNFunc+import Data.Number.ER.RnToRm (FAPWP)+import qualified Data.Number.ER.Real.DomainBox as DBox++import qualified Data.Number.ER.Real as AERN+import Data.Number.ER.Real (IRA)++import Data.Number.ER.Misc++-- function f(x) = x for x in [0,1]:+x :: FAPWP+x =+ AERNFunc.setMaxDegree 2 $+ AERNFunc.proj (DBox.fromAscList [(0,(0) AERN.\/ 1)]) 0+-- function f(x1) = x1 for x1 in [0,1]:+x1 :: FAPWP+x1 =+ AERNFunc.setMaxDegree 2 $+ AERNFunc.proj (DBox.fromAscList [(1,(0) AERN.\/ 1)]) 1++-- domains combined automatically:+fn1 :: FAPWP+fn1 = 2*x + x1++-- ensure the piecewise representation has 4 segments:+fn1depth2 :: FAPWP+fn1depth2 = AERNFunc.bisectUnbisectDepth 2 fn1++-- apply sine pointwise to the function enclosure:+fn2 :: FAPWP+fn2 = +-- AERN.sin 10 fn1depth2+ AERN.sin 15 fn1depth2++-- evaluate the function at point x = 0.1, x1 = 0.1:+fn2at0101 :: IRA+[fn2at0101] = + AERNFunc.eval (DBox.fromList [(0,0.1), (1,0.1)]) fn2++-- partially evaluate fn2 at x1 = 1:+fn3 :: FAPWP+fn3 = AERNFunc.partialEval (DBox.fromList [(1,1)]) fn2++-- integrate fn3 by x with value 1 at origin x = 1:+fn4 :: FAPWP+fn4 = + AERNFunc.integrate ix fn2 var span origin value+ where+ ix = 2 -- effort index+ var = 0+ span = DBox.noinfo -- integrate over the whole domain+ origin = 1+ value = 1++-- integrate fn2 by x1 with value (1 - x) at origin x1 = 0:+fn5 :: FAPWP+fn5 =+ AERNFunc.integrate ix fn2 var span origin value+ where+ ix = 2 -- effort index+ var = 1+ span = DBox.noinfo -- integrate over the whole domain+ origin = 0+ value = 1 - x+++main = + do+ AERN.initMachineDouble+ putStrLn "****************************************"+ putStrLn "Testing polynomial enclosure arithmetic:"+ putStrLn "****************************************"+ putStrLn "**** Projections:"+ putStrLn $+ "x =\n " ++ show x+ putStrLn $+ "\nx1 =\n " ++ show x1+ putStrLn "\n**** Merging domains:"+ putStrLn $+ "2*x + x1 =\n " ++ showHead 12 fn1+ putStrLn "\n**** Bisection depth 2:"+ putStrLn $+ "2*x + x1 =\n " ++ showHead 17 fn1depth2+ putStrLn "\n**** Elementary functions:"+ putStrLn $+ "sin(2*x + x1) =\n " ++ showHead 17 fn2+ putStrLn "\n**** Evaluation:"+ putStrLn $+ "sin(2*x + x1)[x = 0.1, x1 = 0.1] = sin(0.3) = \n " ++ show fn2at0101+ putStrLn "\n**** Partial evaluation:"+ putStrLn $+ "sin(2*x + x1)[x1 = 1] = sin(5*x + 1) = \n " ++ showHead 15 fn3+ putStrLn "\n**** Integration of 1-dim function:"+ putStrLn $+ "f(x) = (Int sin(2*x + 1) dx) [f(1) = 1] =\n " ++ showHead 15 fn4+ putStrLn "\n**** Integration of 2-dim function:"+ putStrLn $+ "f(x,x1) = (Int sin(2*x + x1) dx1) [f(x,1) = 1 - x] =\n " ++ showHead 17 fn5++showHead n = showFirstLastLines n 0