diff --git a/AERN-RnToRm.cabal b/AERN-RnToRm.cabal new file mode 100644 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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