{- |
This module demonstrates triangular matrices.

It verifies that the divided difference scheme
nicely fits into a triangular matrix,
where function addition is mapped to matrix addition
and function multiplication is mapped to matrix multiplication.

<http://en.wikipedia.org/wiki/Divided_difference>
-}
module Numeric.LAPACK.Example.DividedDifference where

import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
import qualified Numeric.LAPACK.Matrix.Square as Square
import qualified Numeric.LAPACK.Matrix.Layout as Layout
import qualified Numeric.LAPACK.Matrix as Matrix
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Matrix (ShapeInt, (#+#), (#-#))
import Numeric.LAPACK.Vector (Vector, (|+|), (|-|))
import Numeric.LAPACK.Format ((##))

import qualified Data.Array.Comfort.Shape as Shape
import qualified Data.Array.Comfort.Storable as Array
import Foreign.Storable (Storable)

import qualified Data.List as List
import Data.Semigroup ((<>))


{- $setup
>>> import qualified Test.Utility as Util
>>> import Test.Utility (approxArray)
>>>
>>> import qualified Numeric.LAPACK.Example.DividedDifference as DD
>>> import qualified Numeric.LAPACK.Vector as Vector
>>> import Numeric.LAPACK.Example.DividedDifference (dividedDifferencesMatrix)
>>> import Numeric.LAPACK.Matrix (ShapeInt, shapeInt, (#+#))
>>> import Numeric.LAPACK.Vector ((|+|))
>>>
>>> import qualified Data.Array.Comfort.Storable as Array
>>>
>>> import qualified Test.QuickCheck as QC
>>>
>>> import Control.Monad (liftM2)
>>> import Data.Tuple.HT (mapPair)
>>> import Data.Semigroup ((<>))
>>>
>>> type Vector = Vector.Vector ShapeInt Float
>>>
>>> genDD :: QC.Gen (Vector, (Vector, Vector))
>>> genDD = do
>>>    (ys0,ys1) <-
>>>       fmap (mapPair (Vector.autoFromList, Vector.autoFromList) .
>>>             unzip . take 10) $
>>>       QC.listOf $ liftM2 (,) (Util.genElement 10) (Util.genElement 10)
>>>    xs <- Util.genDistinct [-10..10] [-10..10] $ Array.shape ys0
>>>    return (xs,(ys0,ys1))
-}


size :: Vector ShapeInt a -> Int
size :: Vector ShapeInt a -> Int
size = ShapeInt -> Int
forall n. ZeroBased n -> n
Shape.zeroBasedSize (ShapeInt -> Int)
-> (Vector ShapeInt a -> ShapeInt) -> Vector ShapeInt a -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector ShapeInt a -> ShapeInt
forall sh a. Array sh a -> sh
Array.shape

subSlices :: Int -> Vector ShapeInt Float -> Vector ShapeInt Float
subSlices :: Int -> Vector ShapeInt Float -> Vector ShapeInt Float
subSlices Int
k Vector ShapeInt Float
xs = Int -> Vector ShapeInt Float -> Vector ShapeInt Float
forall n a.
(Integral n, Storable a) =>
n -> Array (ZeroBased n) a -> Array (ZeroBased n) a
Vector.drop Int
k Vector ShapeInt Float
xs Vector ShapeInt Float
-> Vector ShapeInt Float -> Vector ShapeInt Float
forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
|-| Int -> Vector ShapeInt Float -> Vector ShapeInt Float
forall n a.
(Integral n, Storable a) =>
n -> Array (ZeroBased n) a -> Array (ZeroBased n) a
Vector.take (Vector ShapeInt Float -> Int
forall a. Vector ShapeInt a -> Int
size Vector ShapeInt Float
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k) Vector ShapeInt Float
xs

parameterDifferences :: Vector ShapeInt Float -> [Vector ShapeInt Float]
parameterDifferences :: Vector ShapeInt Float -> [Vector ShapeInt Float]
parameterDifferences Vector ShapeInt Float
xs = (Int -> Vector ShapeInt Float) -> [Int] -> [Vector ShapeInt Float]
forall a b. (a -> b) -> [a] -> [b]
map ((Int -> Vector ShapeInt Float -> Vector ShapeInt Float)
-> Vector ShapeInt Float -> Int -> Vector ShapeInt Float
forall a b c. (a -> b -> c) -> b -> a -> c
flip Int -> Vector ShapeInt Float -> Vector ShapeInt Float
subSlices Vector ShapeInt Float
xs) [Int
1 .. Vector ShapeInt Float -> Int
forall a. Vector ShapeInt a -> Int
size Vector ShapeInt Float
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]

dividedDifferences ::
   Vector ShapeInt Float -> Vector ShapeInt Float -> [Vector ShapeInt Float]
dividedDifferences :: Vector ShapeInt Float
-> Vector ShapeInt Float -> [Vector ShapeInt Float]
dividedDifferences Vector ShapeInt Float
xs Vector ShapeInt Float
ys =
   (Vector ShapeInt Float
 -> Vector ShapeInt Float -> Vector ShapeInt Float)
-> Vector ShapeInt Float
-> [Vector ShapeInt Float]
-> [Vector ShapeInt Float]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl
      (\Vector ShapeInt Float
ddys Vector ShapeInt Float
dxs -> Vector ShapeInt Float
-> Vector ShapeInt Float -> Vector ShapeInt Float
forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
Vector.divide (Int -> Vector ShapeInt Float -> Vector ShapeInt Float
subSlices Int
1 Vector ShapeInt Float
ddys) Vector ShapeInt Float
dxs)
      Vector ShapeInt Float
ys
      (Vector ShapeInt Float -> [Vector ShapeInt Float]
parameterDifferences Vector ShapeInt Float
xs)

upperFromPyramid ::
   (Shape.C sh, Storable a) => sh -> [Vector sh a] -> Triangular.Upper sh a
upperFromPyramid :: sh -> [Vector sh a] -> Upper sh a
upperFromPyramid sh
sh =
   Order -> sh -> [a] -> Upper sh a
forall sh a. (C sh, Storable a) => Order -> sh -> [a] -> Upper sh a
Triangular.upperFromList Order
Layout.RowMajor sh
sh ([a] -> Upper sh a)
-> ([Vector sh a] -> [a]) -> [Vector sh a] -> Upper sh a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[a]] -> [a]) -> ([Vector sh a] -> [[a]]) -> [Vector sh a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[a]] -> [[a]]
forall a. [[a]] -> [[a]]
List.transpose ([[a]] -> [[a]])
-> ([Vector sh a] -> [[a]]) -> [Vector sh a] -> [[a]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector sh a -> [a]) -> [Vector sh a] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map Vector sh a -> [a]
forall sh a. (C sh, Storable a) => Vector sh a -> [a]
Vector.toList

{- |
prop> QC.forAll genDD $ \(xs, (ys0,ys1)) -> approxArray (dividedDifferencesMatrix xs (ys0|+|ys1)) (dividedDifferencesMatrix xs ys0 #+# dividedDifferencesMatrix xs ys1)
prop> QC.forAll genDD $ \(xs, (ys0,ys1)) -> approxArray (dividedDifferencesMatrix xs (Vector.mul ys0 ys1)) (dividedDifferencesMatrix xs ys0 <> dividedDifferencesMatrix xs ys1)
-}
dividedDifferencesMatrix ::
   Vector ShapeInt Float -> Vector ShapeInt Float ->
   Triangular.Upper ShapeInt Float
dividedDifferencesMatrix :: Vector ShapeInt Float
-> Vector ShapeInt Float -> Upper ShapeInt Float
dividedDifferencesMatrix Vector ShapeInt Float
xs Vector ShapeInt Float
ys =
   ShapeInt -> [Vector ShapeInt Float] -> Upper ShapeInt Float
forall sh a.
(C sh, Storable a) =>
sh -> [Vector sh a] -> Upper sh a
upperFromPyramid (Vector ShapeInt Float -> ShapeInt
forall sh a. Array sh a -> sh
Array.shape Vector ShapeInt Float
xs) ([Vector ShapeInt Float] -> Upper ShapeInt Float)
-> [Vector ShapeInt Float] -> Upper ShapeInt Float
forall a b. (a -> b) -> a -> b
$ Vector ShapeInt Float
-> Vector ShapeInt Float -> [Vector ShapeInt Float]
dividedDifferences Vector ShapeInt Float
xs Vector ShapeInt Float
ys


{- |
prop> QC.forAll (QC.choose (0,10)) $ \n -> let sh = shapeInt n in QC.forAll (Util.genDistinct [-10..10] [-10..10] sh) $ \xs -> approxArray (DD.parameterDifferencesMatrix xs) (DD.upperFromPyramid sh (Vector.zero sh : DD.parameterDifferences xs))
-}
parameterDifferencesMatrix ::
   Vector ShapeInt Float -> Triangular.Upper ShapeInt Float
parameterDifferencesMatrix :: Vector ShapeInt Float -> Upper ShapeInt Float
parameterDifferencesMatrix Vector ShapeInt Float
xs =
   let ones :: Vector ShapeInt Float
ones = ShapeInt -> Vector ShapeInt Float
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.one (ShapeInt -> Vector ShapeInt Float)
-> ShapeInt -> Vector ShapeInt Float
forall a b. (a -> b) -> a -> b
$ Vector ShapeInt Float -> ShapeInt
forall sh a. Array sh a -> sh
Array.shape Vector ShapeInt Float
xs
       tp :: Vector ShapeInt Float
-> Vector ShapeInt Float -> General ShapeInt ShapeInt Float
tp = Order
-> Vector ShapeInt Float
-> Vector ShapeInt Float
-> General ShapeInt ShapeInt Float
forall height width a.
(C height, Eq height, C width, Eq width, Floating a) =>
Order
-> Vector height a -> Vector width a -> General height width a
Matrix.tensorProduct Order
Layout.RowMajor
   in Unpacked
  Arbitrary Filled Filled Shape Small Small ShapeInt ShapeInt Float
-> Upper ShapeInt Float
forall property lower meas vert height width a.
(Property property, Strip lower, Measure meas, C vert, C height,
 C width, Floating a) =>
Unpacked property lower Filled meas vert Small height width a
-> Upper width a
Triangular.takeUpper (Unpacked
   Arbitrary Filled Filled Shape Small Small ShapeInt ShapeInt Float
 -> Upper ShapeInt Float)
-> Unpacked
     Arbitrary Filled Filled Shape Small Small ShapeInt ShapeInt Float
-> Upper ShapeInt Float
forall a b. (a -> b) -> a -> b
$ General ShapeInt ShapeInt Float
-> Unpacked
     Arbitrary Filled Filled Shape Small Small ShapeInt ShapeInt Float
forall meas vert horiz sh a.
(Measure meas, C vert, C horiz, Eq sh) =>
Full meas vert horiz sh sh a -> Square sh a
Square.fromFull (General ShapeInt ShapeInt Float
 -> Unpacked
      Arbitrary Filled Filled Shape Small Small ShapeInt ShapeInt Float)
-> General ShapeInt ShapeInt Float
-> Unpacked
     Arbitrary Filled Filled Shape Small Small ShapeInt ShapeInt Float
forall a b. (a -> b) -> a -> b
$ Vector ShapeInt Float
-> Vector ShapeInt Float -> General ShapeInt ShapeInt Float
tp Vector ShapeInt Float
ones Vector ShapeInt Float
xs General ShapeInt ShapeInt Float
-> General ShapeInt ShapeInt Float
-> General ShapeInt ShapeInt Float
forall meas vert horiz typ xl xu height width a lower upper.
(Measure meas, C vert, C horiz, Subtractive typ,
 SubtractiveExtra typ xl, SubtractiveExtra typ xu, C height,
 Eq height, C width, Eq width, Floating a) =>
Matrix typ xl xu lower upper meas vert horiz height width a
-> Matrix typ xl xu lower upper meas vert horiz height width a
-> Matrix typ xl xu lower upper meas vert horiz height width a
#-# Vector ShapeInt Float
-> Vector ShapeInt Float -> General ShapeInt ShapeInt Float
tp Vector ShapeInt Float
xs Vector ShapeInt Float
ones


main :: IO ()
main :: IO ()
main = do
   let xs :: Vector ShapeInt Float
xs  = [Float] -> Vector ShapeInt Float
forall a. Storable a => [a] -> Vector ShapeInt a
Vector.autoFromList [Float
0,Float
1,Float
4,Float
9,Float
16]
   let ys0 :: Vector ShapeInt Float
ys0 = [Float] -> Vector ShapeInt Float
forall a. Storable a => [a] -> Vector ShapeInt a
Vector.autoFromList [Float
3,Float
1,Float
4,Float
1,Float
5]
   let ys1 :: Vector ShapeInt Float
ys1 = [Float] -> Vector ShapeInt Float
forall a. Storable a => [a] -> Vector ShapeInt a
Vector.autoFromList [Float
2,Float
7,Float
1,Float
8,Float
1]

   (Vector ShapeInt Float -> IO ())
-> [Vector ShapeInt Float] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (Vector ShapeInt Float -> String -> IO ()
forall a. Format a => a -> String -> IO ()
## String
"%.4f") ([Vector ShapeInt Float] -> IO ())
-> [Vector ShapeInt Float] -> IO ()
forall a b. (a -> b) -> a -> b
$ Vector ShapeInt Float -> [Vector ShapeInt Float]
parameterDifferences Vector ShapeInt Float
xs
   Vector ShapeInt Float -> Upper ShapeInt Float
parameterDifferencesMatrix Vector ShapeInt Float
xs Upper ShapeInt Float -> String -> IO ()
forall a. Format a => a -> String -> IO ()
## String
"%.4f"

   let ddys0 :: Upper ShapeInt Float
ddys0 = Vector ShapeInt Float
-> Vector ShapeInt Float -> Upper ShapeInt Float
dividedDifferencesMatrix Vector ShapeInt Float
xs Vector ShapeInt Float
ys0
   let ddys1 :: Upper ShapeInt Float
ddys1 = Vector ShapeInt Float
-> Vector ShapeInt Float -> Upper ShapeInt Float
dividedDifferencesMatrix Vector ShapeInt Float
xs Vector ShapeInt Float
ys1
   Upper ShapeInt Float
ddys0 Upper ShapeInt Float -> String -> IO ()
forall a. Format a => a -> String -> IO ()
## String
"%.4f"
   Upper ShapeInt Float
ddys1 Upper ShapeInt Float -> String -> IO ()
forall a. Format a => a -> String -> IO ()
## String
"%.4f"
   String -> IO ()
putStrLn String
""

   Vector ShapeInt Float
-> Vector ShapeInt Float -> Upper ShapeInt Float
dividedDifferencesMatrix Vector ShapeInt Float
xs (Vector ShapeInt Float
ys0Vector ShapeInt Float
-> Vector ShapeInt Float -> Vector ShapeInt Float
forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
|+|Vector ShapeInt Float
ys1) Upper ShapeInt Float -> String -> IO ()
forall a. Format a => a -> String -> IO ()
## String
"%.4f"
   Upper ShapeInt Float
ddys0 Upper ShapeInt Float
-> Upper ShapeInt Float -> Upper ShapeInt Float
forall meas vert horiz typ xl xu height width a lower upper.
(Measure meas, C vert, C horiz, Additive typ, AdditiveExtra typ xl,
 AdditiveExtra typ xu, C height, Eq height, C width, Eq width,
 Floating a) =>
Matrix typ xl xu lower upper meas vert horiz height width a
-> Matrix typ xl xu lower upper meas vert horiz height width a
-> Matrix typ xl xu lower upper meas vert horiz height width a
#+# Upper ShapeInt Float
ddys1 Upper ShapeInt Float -> String -> IO ()
forall a. Format a => a -> String -> IO ()
## String
"%.4f"
   String -> IO ()
putStrLn String
""

   Vector ShapeInt Float
-> Vector ShapeInt Float -> Upper ShapeInt Float
dividedDifferencesMatrix Vector ShapeInt Float
xs (Vector ShapeInt Float
-> Vector ShapeInt Float -> Vector ShapeInt Float
forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
Vector.mul Vector ShapeInt Float
ys0 Vector ShapeInt Float
ys1) Upper ShapeInt Float -> String -> IO ()
forall a. Format a => a -> String -> IO ()
## String
"%.4f"
   Upper ShapeInt Float
ddys0 Upper ShapeInt Float
-> Upper ShapeInt Float -> Upper ShapeInt Float
forall a. Semigroup a => a -> a -> a
<> Upper ShapeInt Float
ddys1 Upper ShapeInt Float -> String -> IO ()
forall a. Format a => a -> String -> IO ()
## String
"%.4f"