{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
-----------------------------------------------------------------------------
-- |
-- Copyright   :  (c) Edward Kmett 2015-2021
-- License     :  BSD3
-- Maintainer  :  ekmett@gmail.com
-- Stability   :  experimental
-- Portability :  GHC only
--
-----------------------------------------------------------------------------

module Numeric.AD.Rank1.Newton.Double
  (
  -- * Newton's Method (Forward)
    findZero
  , inverse
  , fixedPoint
  , extremum
  ) where

import Prelude hiding (all, mapM)
import Numeric.AD.Mode
import Numeric.AD.Rank1.Forward (Forward)
import qualified Numeric.AD.Rank1.Forward as Forward
import Numeric.AD.Rank1.Forward.Double (ForwardDouble, diff')
import Numeric.AD.Internal.On
import Numeric.AD.Internal.Combinators (takeWhileDifferent)

-- | The 'findZero' function finds a zero of a scalar function using
-- Newton's method; its output is a stream of increasingly accurate
-- results.  (Modulo the usual caveats.) If the stream becomes constant
-- ("it converges"), no further elements are returned.
--
-- Examples:
--
-- >>> take 10 $ findZero (\x->x^2-4) 1
-- [1.0,2.5,2.05,2.000609756097561,2.0000000929222947,2.000000000000002,2.0]
findZero :: (ForwardDouble -> ForwardDouble) -> Double -> [Double]
findZero :: (ForwardDouble -> ForwardDouble) -> Double -> [Double]
findZero ForwardDouble -> ForwardDouble
f = [Double] -> [Double]
forall a. Eq a => [a] -> [a]
takeWhileDifferent ([Double] -> [Double])
-> (Double -> [Double]) -> Double -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (ForwardDouble -> ForwardDouble) -> Double -> [Double]
findZeroNoEq ForwardDouble -> ForwardDouble
f
{-# INLINE findZero #-}

-- | The 'findZeroNoEq' function behaves the same as 'findZero' except that it
-- doesn't truncate the list once the results become constant.
findZeroNoEq :: (ForwardDouble -> ForwardDouble) -> Double -> [Double]
findZeroNoEq :: (ForwardDouble -> ForwardDouble) -> Double -> [Double]
findZeroNoEq ForwardDouble -> ForwardDouble
f = [Double] -> [Double]
forall a. Eq a => [a] -> [a]
takeWhileDifferent ([Double] -> [Double])
-> (Double -> [Double]) -> Double -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double) -> Double -> [Double]
forall a. (a -> a) -> a -> [a]
iterate Double -> Double
go where
  go :: Double -> Double
go Double
x = Double
xn where
    (Double
y,Double
y') = (ForwardDouble -> ForwardDouble) -> Double -> (Double, Double)
diff' ForwardDouble -> ForwardDouble
f Double
x
    xn :: Double
xn = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
yDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
y'
{-# INLINE findZeroNoEq #-}

-- | The 'inverse' function inverts a scalar function using
-- Newton's method; its output is a stream of increasingly accurate
-- results.  (Modulo the usual caveats.) If the stream becomes
-- constant ("it converges"), no further elements are returned.
--
-- Example:
--
-- >>> last $ take 10 $ inverse sqrt 1 (sqrt 10)
-- 10.0
inverse :: (ForwardDouble -> ForwardDouble) -> Double -> Double -> [Double]
inverse :: (ForwardDouble -> ForwardDouble) -> Double -> Double -> [Double]
inverse ForwardDouble -> ForwardDouble
f Double
x0 = [Double] -> [Double]
forall a. Eq a => [a] -> [a]
takeWhileDifferent ([Double] -> [Double])
-> (Double -> [Double]) -> Double -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (ForwardDouble -> ForwardDouble) -> Double -> Double -> [Double]
inverseNoEq ForwardDouble -> ForwardDouble
f Double
x0
{-# INLINE inverse  #-}

-- | The 'inverseNoEq' function behaves the same as 'inverse' except that it
-- doesn't truncate the list once the results become constant.
inverseNoEq :: (ForwardDouble -> ForwardDouble) -> Double -> Double -> [Double]
inverseNoEq :: (ForwardDouble -> ForwardDouble) -> Double -> Double -> [Double]
inverseNoEq ForwardDouble -> ForwardDouble
f Double
x0 Double
y = (ForwardDouble -> ForwardDouble) -> Double -> [Double]
findZeroNoEq (\ForwardDouble
x -> ForwardDouble -> ForwardDouble
f ForwardDouble
x ForwardDouble -> ForwardDouble -> ForwardDouble
forall a. Num a => a -> a -> a
- Scalar ForwardDouble -> ForwardDouble
forall t. Mode t => Scalar t -> t
auto Double
Scalar ForwardDouble
y) Double
x0
{-# INLINE inverseNoEq #-}

-- | The 'fixedPoint' function find a fixedpoint of a scalar
-- function using Newton's method; its output is a stream of
-- increasingly accurate results.  (Modulo the usual caveats.)
--
-- If the stream becomes constant ("it converges"), no further
-- elements are returned.
--
-- >>> last $ take 10 $ fixedPoint cos 1
-- 0.7390851332151607
fixedPoint :: (ForwardDouble -> ForwardDouble) -> Double -> [Double]
fixedPoint :: (ForwardDouble -> ForwardDouble) -> Double -> [Double]
fixedPoint ForwardDouble -> ForwardDouble
f = [Double] -> [Double]
forall a. Eq a => [a] -> [a]
takeWhileDifferent ([Double] -> [Double])
-> (Double -> [Double]) -> Double -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (ForwardDouble -> ForwardDouble) -> Double -> [Double]
fixedPointNoEq ForwardDouble -> ForwardDouble
f
{-# INLINE fixedPoint #-}

-- | The 'fixedPointNoEq' function behaves the same as 'fixedPoint' except that
-- doesn't truncate the list once the results become constant.
fixedPointNoEq :: (ForwardDouble -> ForwardDouble) -> Double -> [Double]
fixedPointNoEq :: (ForwardDouble -> ForwardDouble) -> Double -> [Double]
fixedPointNoEq ForwardDouble -> ForwardDouble
f = (ForwardDouble -> ForwardDouble) -> Double -> [Double]
findZeroNoEq (\ForwardDouble
x -> ForwardDouble -> ForwardDouble
f ForwardDouble
x ForwardDouble -> ForwardDouble -> ForwardDouble
forall a. Num a => a -> a -> a
- ForwardDouble
x)
{-# INLINE fixedPointNoEq #-}

-- | The 'extremum' function finds an extremum of a scalar
-- function using Newton's method; produces a stream of increasingly
-- accurate results.  (Modulo the usual caveats.) If the stream
-- becomes constant ("it converges"), no further elements are returned.
--
-- >>> last $ take 10 $ extremum cos 1
-- 0.0
extremum :: (On (Forward ForwardDouble) -> On (Forward ForwardDouble)) -> Double -> [Double]
extremum :: (On (Forward ForwardDouble) -> On (Forward ForwardDouble))
-> Double -> [Double]
extremum On (Forward ForwardDouble) -> On (Forward ForwardDouble)
f = [Double] -> [Double]
forall a. Eq a => [a] -> [a]
takeWhileDifferent ([Double] -> [Double])
-> (Double -> [Double]) -> Double -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (On (Forward ForwardDouble) -> On (Forward ForwardDouble))
-> Double -> [Double]
extremumNoEq On (Forward ForwardDouble) -> On (Forward ForwardDouble)
f
{-# INLINE extremum #-}

-- | The 'extremumNoEq' function behaves the same as 'extremum' except that it
-- doesn't truncate the list once the results become constant.
extremumNoEq :: (On (Forward ForwardDouble) -> On (Forward ForwardDouble)) -> Double -> [Double]
extremumNoEq :: (On (Forward ForwardDouble) -> On (Forward ForwardDouble))
-> Double -> [Double]
extremumNoEq On (Forward ForwardDouble) -> On (Forward ForwardDouble)
f = (ForwardDouble -> ForwardDouble) -> Double -> [Double]
findZeroNoEq ((Forward ForwardDouble -> Forward ForwardDouble)
-> ForwardDouble -> ForwardDouble
forall a. Num a => (Forward a -> Forward a) -> a -> a
Forward.diff (On (Forward ForwardDouble) -> Forward ForwardDouble
forall t. On t -> t
off (On (Forward ForwardDouble) -> Forward ForwardDouble)
-> (Forward ForwardDouble -> On (Forward ForwardDouble))
-> Forward ForwardDouble
-> Forward ForwardDouble
forall b c a. (b -> c) -> (a -> b) -> a -> c
. On (Forward ForwardDouble) -> On (Forward ForwardDouble)
f (On (Forward ForwardDouble) -> On (Forward ForwardDouble))
-> (Forward ForwardDouble -> On (Forward ForwardDouble))
-> Forward ForwardDouble
-> On (Forward ForwardDouble)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Forward ForwardDouble -> On (Forward ForwardDouble)
forall t. t -> On t
On))
{-# INLINE extremumNoEq #-}