summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorMichalKonecny <>2008-08-07 10:52:03 (GMT)
committerLuite Stegeman <luite@luite.com>2008-08-07 10:52:03 (GMT)
commitb04c2547140974dd236a008d1a7734c9e673596c (patch)
treec2ec0e4c222f697e26bb702a1cec03022e3d0eab
version 0.3.00.3.0
-rw-r--r--AERN-RnToRm.cabal123
-rw-r--r--ChangeLog4
-rw-r--r--LICENCE30
-rw-r--r--Setup.lhs3
-rw-r--r--src/Data/Number/ER/RnToRm.hs100
-rw-r--r--src/Data/Number/ER/RnToRm/Approx.hs303
-rw-r--r--src/Data/Number/ER/RnToRm/Approx/DomEdges.hs487
-rw-r--r--src/Data/Number/ER/RnToRm/Approx/DomTransl.hs495
-rw-r--r--src/Data/Number/ER/RnToRm/Approx/PieceWise.hs522
-rw-r--r--src/Data/Number/ER/RnToRm/Approx/Tuple.hs272
-rw-r--r--src/Data/Number/ER/RnToRm/BisectionTree.hs665
-rw-r--r--src/Data/Number/ER/RnToRm/BisectionTree/Integration.hs278
-rw-r--r--src/Data/Number/ER/RnToRm/BisectionTree/Path.hs134
-rw-r--r--src/Data/Number/ER/RnToRm/DefaultRepr.hs64
-rw-r--r--src/Data/Number/ER/RnToRm/TestingDefs.hs72
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/Approx.hs92
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs589
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/Base.hs348
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs93
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs279
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs293
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs455
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs190
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs226
-rw-r--r--src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs169
-rw-r--r--tests/Demo.hs114
26 files changed, 6400 insertions, 0 deletions
diff --git a/AERN-RnToRm.cabal b/AERN-RnToRm.cabal
new file mode 100644
index 0000000..3625f7b
--- /dev/null
+++ b/AERN-RnToRm.cabal
@@ -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
+
diff --git a/ChangeLog b/ChangeLog
new file mode 100644
index 0000000..27a3584
--- /dev/null
+++ b/ChangeLog
@@ -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
+
diff --git a/LICENCE b/LICENCE
new file mode 100644
index 0000000..ac611a5
--- /dev/null
+++ b/LICENCE
@@ -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.
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
index 0000000..5bde0de
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#!/usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/src/Data/Number/ER/RnToRm.hs b/src/Data/Number/ER/RnToRm.hs
new file mode 100644
index 0000000..cb87ff7
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm.hs
@@ -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
+
diff --git a/src/Data/Number/ER/RnToRm/Approx.hs b/src/Data/Number/ER/RnToRm/Approx.hs
new file mode 100644
index 0000000..7f41bba
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/Approx.hs
@@ -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
diff --git a/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs b/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs
new file mode 100644
index 0000000..b5176fe
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/Approx/DomEdges.hs
@@ -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
+
+
diff --git a/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs b/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs
new file mode 100644
index 0000000..780770f
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/Approx/DomTransl.hs
@@ -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
+
+
diff --git a/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs b/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs
new file mode 100644
index 0000000..968b506
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/Approx/PieceWise.hs
@@ -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
+
diff --git a/src/Data/Number/ER/RnToRm/Approx/Tuple.hs b/src/Data/Number/ER/RnToRm/Approx/Tuple.hs
new file mode 100644
index 0000000..4b37c15
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/Approx/Tuple.hs
@@ -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
diff --git a/src/Data/Number/ER/RnToRm/BisectionTree.hs b/src/Data/Number/ER/RnToRm/BisectionTree.hs
new file mode 100644
index 0000000..093a854
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/BisectionTree.hs
@@ -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
+
diff --git a/src/Data/Number/ER/RnToRm/BisectionTree/Integration.hs b/src/Data/Number/ER/RnToRm/BisectionTree/Integration.hs
new file mode 100644
index 0000000..f7dc5ba
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/BisectionTree/Integration.hs
@@ -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)
+
diff --git a/src/Data/Number/ER/RnToRm/BisectionTree/Path.hs b/src/Data/Number/ER/RnToRm/BisectionTree/Path.hs
new file mode 100644
index 0000000..9bd6ea7
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/BisectionTree/Path.hs
@@ -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
diff --git a/src/Data/Number/ER/RnToRm/DefaultRepr.hs b/src/Data/Number/ER/RnToRm/DefaultRepr.hs
new file mode 100644
index 0000000..ff33ada
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/DefaultRepr.hs
@@ -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
+
diff --git a/src/Data/Number/ER/RnToRm/TestingDefs.hs b/src/Data/Number/ER/RnToRm/TestingDefs.hs
new file mode 100644
index 0000000..3b1b57f
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/TestingDefs.hs
@@ -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!
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/Approx.hs b/src/Data/Number/ER/RnToRm/UnitDom/Approx.hs
new file mode 100644
index 0000000..1e0b989
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/UnitDom/Approx.hs
@@ -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
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs b/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs
new file mode 100644
index 0000000..86bdc3f
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/UnitDom/Approx/Interval.hs
@@ -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))
+
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/Base.hs b/src/Data/Number/ER/RnToRm/UnitDom/Base.hs
new file mode 100644
index 0000000..0a041b3
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/UnitDom/Base.hs
@@ -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
+
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs
new file mode 100644
index 0000000..809f92b
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom.hs
@@ -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
+
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs
new file mode 100644
index 0000000..4ede2ec
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Basic.hs
@@ -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)
+
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs
new file mode 100644
index 0000000..68b2699
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Bounds.hs
@@ -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
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs
new file mode 100644
index 0000000..20b80fd
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Elementary.hs
@@ -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
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs
new file mode 100644
index 0000000..7907aac
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Eval.hs
@@ -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
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs
new file mode 100644
index 0000000..8fd15f1
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Field.hs
@@ -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"
+
diff --git a/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs
new file mode 100644
index 0000000..461103b
--- /dev/null
+++ b/src/Data/Number/ER/RnToRm/UnitDom/ChebyshevBase/Polynom/Integration.hs
@@ -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"
+
diff --git a/tests/Demo.hs b/tests/Demo.hs
new file mode 100644
index 0000000..97f483b
--- /dev/null
+++ b/tests/Demo.hs
@@ -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