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 ((<>))
size :: Vector ShapeInt a -> Int
size = Shape.zeroBasedSize . Array.shape
subSlices :: Int -> Vector ShapeInt Float -> Vector ShapeInt Float
subSlices k xs = Vector.drop k xs |-| Vector.take (size xs - k) xs
parameterDifferences :: Vector ShapeInt Float -> [Vector ShapeInt Float]
parameterDifferences xs = map (flip subSlices xs) [1 .. size xs - 1]
dividedDifferences ::
Vector ShapeInt Float -> Vector ShapeInt Float -> [Vector ShapeInt Float]
dividedDifferences xs ys =
scanl
(\ddys dxs -> Vector.divide (subSlices 1 ddys) dxs)
ys
(parameterDifferences xs)
upperFromPyramid ::
(Shape.C sh, Storable a) => sh -> [Vector sh a] -> Triangular.Upper sh a
upperFromPyramid sh =
Triangular.upperFromList Layout.RowMajor sh .
concat . List.transpose . map Vector.toList
dividedDifferencesMatrix ::
Vector ShapeInt Float -> Vector ShapeInt Float ->
Triangular.Upper ShapeInt Float
dividedDifferencesMatrix xs ys =
upperFromPyramid (Array.shape xs) $ dividedDifferences xs ys
parameterDifferencesMatrix ::
Vector ShapeInt Float -> Triangular.Upper ShapeInt Float
parameterDifferencesMatrix xs =
let ones = Vector.one $ Array.shape xs
tp = Matrix.tensorProduct Layout.RowMajor
in Triangular.takeUpper $ Square.fromFull $ tp ones xs #-# tp xs ones
main :: IO ()
main = do
let xs = Vector.autoFromList [0,1,4,9,16]
let ys0 = Vector.autoFromList [3,1,4,1,5]
let ys1 = Vector.autoFromList [2,7,1,8,1]
mapM_ (## "%.4f") $ parameterDifferences xs
parameterDifferencesMatrix xs ## "%.4f"
let ddys0 = dividedDifferencesMatrix xs ys0
let ddys1 = dividedDifferencesMatrix xs ys1
ddys0 ## "%.4f"
ddys1 ## "%.4f"
putStrLn ""
dividedDifferencesMatrix xs (ys0|+|ys1) ## "%.4f"
ddys0 #+# ddys1 ## "%.4f"
putStrLn ""
dividedDifferencesMatrix xs (Vector.mul ys0 ys1) ## "%.4f"
ddys0 <> ddys1 ## "%.4f"