{-# LANGUAGE RebindableSyntax #-}
module MathObj.RefinementMask2 (
   T, coeffs, fromCoeffs,
   fromPolynomial,
   toPolynomial,
   toPolynomialFast,
   refinePolynomial,
   convolvePolynomial,
   convolveTruncatedPowerPolynomials,
   ) where

import qualified MathObj.Polynomial as Poly
import qualified MathObj.Polynomial.Core as PolyCore
import qualified Algebra.RealField as RealField
import qualified Algebra.Field  as Field
import qualified Algebra.Ring   as Ring
import qualified Algebra.Vector as Vector

import qualified Data.List as List
import qualified Data.List.HT as ListHT
import qualified Data.List.Match as Match
import Data.Maybe (fromMaybe, )
import Control.Monad (liftM2, )

import qualified Test.QuickCheck as QC

import qualified NumericPrelude.List.Generic as NPList
import NumericPrelude.Base
import NumericPrelude.Numeric


{- $setup
>>> import qualified MathObj.RefinementMask2 as Mask
>>> import qualified MathObj.Polynomial      as Poly
>>> import qualified MathObj.Polynomial.Core as PolyCore
>>>
>>> import qualified Algebra.Differential as D
>>> import qualified Algebra.Ring as Ring
>>> import Test.NumericPrelude.Utility ((/\))
>>> import qualified Test.QuickCheck as QC
>>> import NumericPrelude.Numeric as NP
>>> import NumericPrelude.Base as P
>>> import Prelude ()
>>>
>>> import Data.Function.HT (nest)
>>> import Data.Maybe (fromMaybe)
>>>
>>>
>>> hasMultipleZero :: (Ring.C a, Eq a) => Int -> a -> Poly.T a -> Bool
>>> hasMultipleZero n x poly =
>>>    all (zero==) $ take n $
>>>    map (flip Poly.evaluate x) $
>>>    iterate D.differentiate poly
>>>
>>> genAdmissibleMask :: QC.Gen (Mask.T Rational, Poly.T Rational)
>>> genAdmissibleMask =
>>>    QC.suchThatMap QC.arbitrary $
>>>       \mask -> fmap ((,) mask) $ Mask.toPolynomial mask
>>>
>>> polyFromMask :: Mask.T a -> Poly.T a
>>> polyFromMask = Poly.fromCoeffs . Mask.coeffs
>>>
>>> genShortPolynomial :: Int -> QC.Gen (Poly.T Rational)
>>> genShortPolynomial n =
>>>    fmap (Poly.fromCoeffs . PolyCore.normalize . take n) $ QC.arbitrary
-}


newtype T a = Cons {T a -> [a]
coeffs :: [a]}


{-# INLINE fromCoeffs #-}
fromCoeffs :: [a] -> T a
fromCoeffs :: [a] -> T a
fromCoeffs = [a] -> T a
forall a. [a] -> T a
lift0

{-# INLINE lift0 #-}
lift0 :: [a] -> T a
lift0 :: [a] -> T a
lift0 = [a] -> T a
forall a. [a] -> T a
Cons

{-
{-# INLINE lift1 #-}
lift1 :: ([a] -> [a]) -> (T a -> T a)
lift1 f (Cons x0) = Cons (f x0)

{-# INLINE lift2 #-}
lift2 :: ([a] -> [a] -> [a]) -> (T a -> T a -> T a)
lift2 f (Cons x0) (Cons x1) = Cons (f x0 x1)
-}

{-
Functor instance is e.g. useful for converting number types,
say 'Rational' to 'Double'.
-}

instance Functor T where
   fmap :: (a -> b) -> T a -> T b
fmap a -> b
f (Cons [a]
xs) = [b] -> T b
forall a. [a] -> T a
Cons ((a -> b) -> [a] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map a -> b
f [a]
xs)

{-# INLINE appPrec #-}
appPrec :: Int
appPrec :: Int
appPrec  = Int
10

instance (Show a) => Show (T a) where
   showsPrec :: Int -> T a -> ShowS
showsPrec Int
p (Cons [a]
xs) =
      Bool -> ShowS -> ShowS
showParen (Int
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
appPrec)
         (String -> ShowS
showString String
"RefinementMask2.fromCoeffs " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> ShowS
forall a. Show a => a -> ShowS
shows [a]
xs)

instance (QC.Arbitrary a, Field.C a) => QC.Arbitrary (T a) where
   arbitrary :: Gen (T a)
arbitrary =
      (Integer -> [a] -> T a) -> Gen Integer -> Gen [a] -> Gen (T a)
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2
         (\Integer
degree [a]
body ->
            let s :: a
s = [a] -> a
forall a. C a => [a] -> a
sum [a]
body
            in  [a] -> T a
forall a. [a] -> T a
Cons ([a] -> T a) -> [a] -> T a
forall a b. (a -> b) -> a -> b
$ (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ((a
2 a -> Integer -> a
forall a. C a => a -> Integer -> a
^- Integer
degree a -> a -> a
forall a. C a => a -> a -> a
- a
s) a -> a -> a
forall a. C a => a -> a -> a
/ [a] -> a
forall n a. C n => [a] -> n
NPList.lengthLeft [a]
body a -> a -> a
forall a. C a => a -> a -> a
+) [a]
body)
         ((Integer, Integer) -> Gen Integer
forall a. Random a => (a, a) -> Gen a
QC.choose (-Integer
5,Integer
0)) Gen [a]
forall a. Arbitrary a => Gen a
QC.arbitrary


{- |
Determine mask by Gauss elimination.

R - alternating binomial coefficients
L - differences of translated polynomials in columns

p2 = L * R^(-1) * m

R * L^(-1) * p2 = m


prop> genAdmissibleMask /\ \(mask,poly) -> hasMultipleZero (fromMaybe 0 $ Poly.degree poly) 1 (polyFromMask (Mask.fromPolynomial poly) - polyFromMask mask)

prop> genShortPolynomial 5 /\ \poly -> maybe False (Poly.collinear poly) $ Mask.toPolynomial $ Mask.fromPolynomial poly
-}
fromPolynomial ::
   (Field.C a) => Poly.T a -> T a
fromPolynomial :: T a -> T a
fromPolynomial T a
poly =
   [a] -> T a
forall a. [a] -> T a
fromCoeffs ([a] -> T a) -> [a] -> T a
forall a b. (a -> b) -> a -> b
$
   (a -> [a] -> [a]) -> [a] -> [a] -> [a]
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\a
p [a]
ps ->
      (a -> a -> a) -> [a] -> [a]
forall a b. (a -> a -> b) -> [a] -> [b]
ListHT.mapAdjacent (-) (a
pa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
ps[a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++[a
0]))
      [] ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$
   ((a, [a]) -> [a] -> [a]) -> [a] -> [(a, [a])] -> [a]
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\(a
db,[a]
dp) [a]
cs ->
      [a] -> ([a] -> a -> [a]) -> [a] -> [a]
forall b a. b -> ([a] -> a -> b) -> [a] -> b
ListHT.switchR
         (String -> [a]
forall a. HasCallStack => String -> a
error String
"RefinementMask2.fromPolynomial: polynomial should be non-empty")
         (\[a]
dps a
dpe ->
            [a]
cs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [(a
db a -> a -> a
forall a. C a => a -> a -> a
- [a] -> [a] -> a
forall a. C a => [a] -> [a] -> a
Ring.scalarProduct [a]
dps [a]
cs) a -> a -> a
forall a. C a => a -> a -> a
/ a
dpe])
         [a]
dp) [] ([(a, [a])] -> [a]) -> [(a, [a])] -> [a]
forall a b. (a -> b) -> a -> b
$
   [a] -> [[a]] -> [(a, [a])]
forall a b. [a] -> [b] -> [(a, b)]
zip
      (T a -> [a]
forall a. T a -> [a]
Poly.coeffs (T a -> [a]) -> T a -> [a]
forall a b. (a -> b) -> a -> b
$ a -> T a -> T a
forall a. C a => a -> T a -> T a
Poly.dilate a
2 T a
poly)
      ([[a]] -> [[a]]
forall a. [[a]] -> [[a]]
List.transpose ([[a]] -> [[a]]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> a -> b
$
       [a] -> [[a]] -> [[a]]
forall b a. [b] -> [a] -> [a]
Match.take (T a -> [a]
forall a. T a -> [a]
Poly.coeffs T a
poly) ([[a]] -> [[a]]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> a -> b
$
       (T a -> [a]) -> [T a] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map T a -> [a]
forall a. T a -> [a]
Poly.coeffs ([T a] -> [[a]]) -> [T a] -> [[a]]
forall a b. (a -> b) -> a -> b
$
       (T a -> T a) -> T a -> [T a]
forall a. (a -> a) -> a -> [a]
iterate T a -> T a
forall a. C a => T a -> T a
polynomialDifference T a
poly)

polynomialDifference ::
   (Ring.C a) => Poly.T a -> Poly.T a
polynomialDifference :: T a -> T a
polynomialDifference T a
poly =
   [a] -> T a
forall a. [a] -> T a
Poly.fromCoeffs ([a] -> T a) -> [a] -> T a
forall a b. (a -> b) -> a -> b
$ [a] -> [a]
forall a. [a] -> [a]
init ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ T a -> [a]
forall a. T a -> [a]
Poly.coeffs (T a -> [a]) -> T a -> [a]
forall a b. (a -> b) -> a -> b
$
   a -> T a -> T a
forall a. C a => a -> T a -> T a
Poly.translate a
1 T a
poly T a -> T a -> T a
forall a. C a => a -> a -> a
- T a
poly

{- |
If the mask does not sum up to a power of @1/2@
then the function returns 'Nothing'.

>>> fmap ((6::Rational) *>) $ Mask.toPolynomial (Mask.fromCoeffs [0.1, 0.02, 0.005::Rational])
Just (Polynomial.fromCoeffs [-12732 % 109375,272 % 625,-18 % 25,1 % 1])
-}
toPolynomial ::
   (RealField.C a) => T a -> Maybe (Poly.T a)
toPolynomial :: T a -> Maybe (T a)
toPolynomial (Cons []) = T a -> Maybe (T a)
forall a. a -> Maybe a
Just (T a -> Maybe (T a)) -> T a -> Maybe (T a)
forall a b. (a -> b) -> a -> b
$ [a] -> T a
forall a. [a] -> T a
Poly.fromCoeffs []
toPolynomial T a
mask =
   let s :: a
s = [a] -> a
forall a. C a => [a] -> a
sum ([a] -> a) -> [a] -> a
forall a b. (a -> b) -> a -> b
$ T a -> [a]
forall a. T a -> [a]
coeffs T a
mask
       ks :: [a]
ks = [a] -> [a]
forall a. [a] -> [a]
reverse ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<=a
1) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a
2a -> a -> a
forall a. C a => a -> a -> a
*) a
s
   in  case [a]
ks of
          a
1:[a]
ks0 ->
             T a -> Maybe (T a)
forall a. a -> Maybe a
Just (T a -> Maybe (T a)) -> T a -> Maybe (T a)
forall a b. (a -> b) -> a -> b
$
             (T a -> a -> T a) -> T a -> [a] -> T a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                (\T a
p a
k ->
                   let ip :: T a
ip = a -> T a -> T a
forall a. C a => a -> T a -> T a
Poly.integrate a
forall a. C a => a
zero T a
p
                   in  T a
ip T a -> T a -> T a
forall a. C a => a -> a -> a
+ a -> T a
forall a. a -> T a
Poly.const (T a -> T a -> a
forall a. C a => T a -> T a -> a
correctConstant ((a -> a) -> T a -> T a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a
ka -> a -> a
forall a. C a => a -> a -> a
/a
sa -> a -> a
forall a. C a => a -> a -> a
*) T a
mask) T a
ip))
                (a -> T a
forall a. a -> T a
Poly.const a
1) [a]
ks0
          [a]
_ -> Maybe (T a)
forall a. Maybe a
Nothing

{-
The constant term must be zero,
higher terms must already satisfy the refinement constraint.
-}
correctConstant ::
   (Field.C a) => T a -> Poly.T a -> a
correctConstant :: T a -> T a -> a
correctConstant T a
mask T a
poly =
   let refined :: T a
refined = T a -> T a -> T a
forall a. C a => T a -> T a -> T a
refinePolynomial T a
mask T a
poly
   in  [a] -> a
forall a. [a] -> a
head (T a -> [a]
forall a. T a -> [a]
Poly.coeffs T a
refined) a -> a -> a
forall a. C a => a -> a -> a
/ (a
1 a -> a -> a
forall a. C a => a -> a -> a
- [a] -> a
forall a. C a => [a] -> a
sum (T a -> [a]
forall a. T a -> [a]
coeffs T a
mask))

toPolynomialFast ::
   (RealField.C a) => T a -> Maybe (Poly.T a)
toPolynomialFast :: T a -> Maybe (T a)
toPolynomialFast T a
mask =
   let s :: a
s = [a] -> a
forall a. C a => [a] -> a
sum ([a] -> a) -> [a] -> a
forall a b. (a -> b) -> a -> b
$ T a -> [a]
forall a. T a -> [a]
coeffs T a
mask
       ks :: [a]
ks = [a] -> [a]
forall a. [a] -> [a]
reverse ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<=a
1) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a
2a -> a -> a
forall a. C a => a -> a -> a
*) a
s
   in  case [a]
ks of
          a
1:[a]
ks0 ->
             T a -> Maybe (T a)
forall a. a -> Maybe a
Just (T a -> Maybe (T a)) -> T a -> Maybe (T a)
forall a b. (a -> b) -> a -> b
$
             (T a -> a -> T a) -> T a -> [a] -> T a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                (\T a
p a
k ->
                   let ip :: T a
ip = a -> T a -> T a
forall a. C a => a -> T a -> T a
Poly.integrate a
forall a. C a => a
zero T a
p
                       c :: a
c = [a] -> a
forall a. [a] -> a
head (T a -> [a]
forall a. T a -> [a]
Poly.coeffs (T a -> T a -> T a
forall a. C a => T a -> T a -> T a
refinePolynomial T a
mask T a
ip))
                   in  T a
ip T a -> T a -> T a
forall a. C a => a -> a -> a
+ a -> T a
forall a. a -> T a
Poly.const (a
ca -> a -> a
forall a. C a => a -> a -> a
*a
k a -> a -> a
forall a. C a => a -> a -> a
/ ((a
1a -> a -> a
forall a. C a => a -> a -> a
-a
k)a -> a -> a
forall a. C a => a -> a -> a
*a
s)))
                (a -> T a
forall a. a -> T a
Poly.const a
1) [a]
ks0
          [a]
_ -> Maybe (T a)
forall a. Maybe a
Nothing

{- |
prop> genShortPolynomial 5 /\ \poly -> poly == Mask.refinePolynomial (Mask.fromPolynomial poly) poly

>>> fmap (round :: Double -> Integer) $ fmap (1000000*) $ nest 50 (Mask.refinePolynomial (Mask.fromCoeffs [0.1, 0.02, 0.005])) (Poly.fromCoeffs [0,0,0,1])
Polynomial.fromCoeffs [-116407,435200,-720000,1000000]
-}
refinePolynomial ::
   (Ring.C a) => T a -> Poly.T a -> Poly.T a
refinePolynomial :: T a -> T a -> T a
refinePolynomial T a
mask =
   a -> T a -> T a
forall a. C a => a -> T a -> T a
Poly.shrink a
2 (T a -> T a) -> (T a -> T a) -> T a -> T a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   [a] -> [T a] -> T a
forall a (v :: * -> *). (C a, C v) => [a] -> [v a] -> v a
Vector.linearComb (T a -> [a]
forall a. T a -> [a]
coeffs T a
mask) ([T a] -> T a) -> (T a -> [T a]) -> T a -> T a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (T a -> T a) -> T a -> [T a]
forall a. (a -> a) -> a -> [a]
iterate (a -> T a -> T a
forall a. C a => a -> T a -> T a
Poly.translate a
1)

convolve ::
   (Ring.C a) => T a -> T a -> T a
convolve :: T a -> T a -> T a
convolve T a
x T a
y =
   [a] -> T a
forall a. [a] -> T a
fromCoeffs ([a] -> T a) -> [a] -> T a
forall a b. (a -> b) -> a -> b
$
   [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PolyCore.mul (T a -> [a]
forall a. T a -> [a]
coeffs T a
x) (T a -> [a]
forall a. T a -> [a]
coeffs T a
y)

{- |
Convolve polynomials via refinement mask.

(mask x + ux*(-1,1)^degree x) * (mask y + uy*(-1,1)^degree y)
-}
convolvePolynomial ::
   (RealField.C a) =>
   Poly.T a -> Poly.T a -> Poly.T a
convolvePolynomial :: T a -> T a -> T a
convolvePolynomial T a
x T a
y =
   T a -> Maybe (T a) -> T a
forall a. a -> Maybe a -> a
fromMaybe
      (String -> T a
forall a. HasCallStack => String -> a
error String
"RefinementMask2.convolvePolynomial: leading term should always be correct") (Maybe (T a) -> T a) -> Maybe (T a) -> T a
forall a b. (a -> b) -> a -> b
$
   T a -> Maybe (T a)
forall a. C a => T a -> Maybe (T a)
toPolynomial (T a -> Maybe (T a)) -> T a -> Maybe (T a)
forall a b. (a -> b) -> a -> b
$ (a -> a) -> T a -> T a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a -> a -> a
forall a. C a => a -> a -> a
/a
2) (T a -> T a) -> T a -> T a
forall a b. (a -> b) -> a -> b
$
   T a -> T a -> T a
forall a. C a => T a -> T a -> T a
convolve (T a -> T a
forall a. C a => T a -> T a
fromPolynomial T a
x) (T a -> T a
forall a. C a => T a -> T a
fromPolynomial T a
y)

{-
This function interprets all monomials as truncated power functions,
that is power functions that are set to zero for negative arguments.
However the convolution implied by this interpretation
cannot be represented by means of mask convolution.
See for instance:

*MathObj.RefinementMask2> let x = Poly.fromCoeffs [1,1] :: Poly.T Rational
*MathObj.RefinementMask2> fromPolynomial $ convolvePolynomial2 x x
RefinementMask2.fromCoeffs [1 % 3,-1 % 8,-1 % 8,1 % 24]

The obtained mask cannot be factored,
thus it is not a complete square.
But maybe it becomes a square if we add u*(-1,1)^4.
However this mask has sum 1/8 and the added term has sum 0,
thus the sum of the modified mask is still 1/8 and thus not a square.
-}
convolveTruncatedPowerPolynomials ::
   (RealField.C a) =>
   Poly.T a -> Poly.T a -> Poly.T a
convolveTruncatedPowerPolynomials :: T a -> T a -> T a
convolveTruncatedPowerPolynomials T a
x T a
y =
   let facs :: [a]
facs = (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl a -> a -> a
forall a. C a => a -> a -> a
(*) a
1 ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a
1a -> a -> a
forall a. C a => a -> a -> a
+) a
1
       xl :: [a]
xl = T a -> [a]
forall a. T a -> [a]
Poly.coeffs T a
x
       yl :: [a]
yl = T a -> [a]
forall a. T a -> [a]
Poly.coeffs T a
y
   in  a -> T a -> T a
forall a. C a => a -> T a -> T a
Poly.integrate a
0 (T a -> T a) -> T a -> T a
forall a b. (a -> b) -> a -> b
$
       [a] -> T a
forall a. [a] -> T a
Poly.fromCoeffs ([a] -> T a) -> [a] -> T a
forall a b. (a -> b) -> a -> b
$
       (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((a -> a -> a) -> a -> a -> a
forall a b c. (a -> b -> c) -> b -> a -> c
flip a -> a -> a
forall a. C a => a -> a -> a
(/)) [a]
facs ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$
       [a] -> [a] -> [a]
forall a. C a => [a] -> [a] -> [a]
PolyCore.mul
          ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. C a => a -> a -> a
(*) [a]
facs [a]
xl)
          ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. C a => a -> a -> a
(*) [a]
facs [a]
yl)