{- |
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.Shape as MatrixShape
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 qualified Data.List as List
import Data.Semigroup ((<>))


{- $setup
>>> import qualified Test.Utility as Util
>>> import Test.Utility (approxArray)
>>>
>>> import qualified Numeric.LAPACK.Vector as Vector
>>> import Numeric.LAPACK.Example.DividedDifference (dividedDifferencesMatrix)
>>> import Numeric.LAPACK.Matrix (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 $ 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)

{- |
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 =
   Order -> ShapeInt -> [Float] -> Upper ShapeInt Float
forall sh a. (C sh, Storable a) => Order -> sh -> [a] -> Upper sh a
Triangular.upperFromList Order
MatrixShape.RowMajor (Vector ShapeInt Float -> ShapeInt
forall sh a. Array sh a -> sh
Array.shape Vector ShapeInt Float
xs) ([Float] -> Upper ShapeInt Float)
-> [Float] -> Upper ShapeInt Float
forall a b. (a -> b) -> a -> b
$
   [[Float]] -> [Float]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[Float]] -> [Float]) -> [[Float]] -> [Float]
forall a b. (a -> b) -> a -> b
$ [[Float]] -> [[Float]]
forall a. [[a]] -> [[a]]
List.transpose ([[Float]] -> [[Float]]) -> [[Float]] -> [[Float]]
forall a b. (a -> b) -> a -> b
$ (Vector ShapeInt Float -> [Float])
-> [Vector ShapeInt Float] -> [[Float]]
forall a b. (a -> b) -> [a] -> [b]
map Vector ShapeInt Float -> [Float]
forall sh a. (C sh, Storable a) => Vector sh a -> [a]
Vector.toList ([Vector ShapeInt Float] -> [[Float]])
-> [Vector ShapeInt Float] -> [[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


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
MatrixShape.RowMajor
   in Full Small Small ShapeInt ShapeInt Float -> Upper ShapeInt Float
forall vert height width a.
(C vert, C height, C width, Floating a) =>
Full vert Small height width a -> Upper width a
Triangular.takeUpper (Full Small Small ShapeInt ShapeInt Float -> Upper ShapeInt Float)
-> Full Small Small ShapeInt ShapeInt Float -> Upper ShapeInt Float
forall a b. (a -> b) -> a -> b
$ General ShapeInt ShapeInt Float
-> Full Small Small ShapeInt ShapeInt Float
forall sh a. Eq sh => General sh sh a -> Square sh a
Square.fromGeneral (General ShapeInt ShapeInt Float
 -> Full Small Small ShapeInt ShapeInt Float)
-> General ShapeInt ShapeInt Float
-> Full 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 shape a.
(Additive shape, Floating a) =>
ArrayMatrix shape a -> ArrayMatrix shape a -> ArrayMatrix shape 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 shape a.
(Additive shape, Floating a) =>
ArrayMatrix shape a -> ArrayMatrix shape a -> ArrayMatrix shape 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"