Basic.hs

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