ad 4.2 → 4.2.0.1

raw patch · 3 files changed

+105/−11 lines, 3 filesdep ~base

Dependency ranges changed: base

Files

• ad.cabal +4/−3
• bench/BlackScholes.hs +12/−8
• src/Numeric/AD/Rank1/Halley.hs +89/−0

@@ -1,5 +1,5 @@ name:         ad-version:      4.2+version:      4.2.0.1 license:      BSD3 license-File: LICENSE copyright:    (c) Edward Kmett 2010-2014,@@ -142,10 +142,11 @@     Numeric.AD.Newton     Numeric.AD.Rank1.Forward     Numeric.AD.Rank1.Forward.Double-    Numeric.AD.Rank1.Sparse-    Numeric.AD.Rank1.Tower+    Numeric.AD.Rank1.Halley     Numeric.AD.Rank1.Kahn     Numeric.AD.Rank1.Newton+    Numeric.AD.Rank1.Sparse+    Numeric.AD.Rank1.Tower    other-modules:     Numeric.AD.Internal.Combinators

@@ -3,6 +3,7 @@ import Data.Number.Erf import qualified Numeric.AD as Mixed import qualified Numeric.AD.Mode.Forward as Forward+import qualified Numeric.AD.Mode.Forward.Double as ForwardDouble import qualified Numeric.AD.Mode.Kahn as Kahn import qualified Numeric.AD.Mode.Reverse as Reverse import qualified Numeric.AD.Mode.Sparse as Sparse@@ -39,6 +40,9 @@         [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Forward.jacobian (fromPair . bs) [r, s, v, t, k]) 2         , bench "greeks Float" $ nf (runFloat $ \r s v t k -> Forward.jacobian (fromPair . bs) [r, s, v, t, k]) 2         ]+    , bgroup "ForwardDouble"+        [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> ForwardDouble.jacobian (fromPair . bs) [r, s, v, t, k]) 2+        ]     , bgroup "Kahn"         [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Kahn.jacobian (fromPair . bs) [r, s, v, t, k]) 2         , bench "higherGreeks Double" $ nf (runDouble $ \r s v t k -> (Kahn.hessian (fst . bs) [r, s, v, t, k], Kahn.hessian (snd . bs) [r, s, v, t, k])) 2@@ -63,12 +67,12 @@         , bench "higherGreeks Float" $ nf (runFloat $ \r s v t k -> (Sparse.hessian (fst . bs) [r, s, v, t, k], Sparse.hessian (snd . bs) [r, s, v, t, k])) 2         , bench "highererGreeks Float" $ nf (runFloat $ \r s v t k -> Sparse.hessianF (fromPair . bs) [r, s, v, t, k]) 2         ]---    , bgroup "Mixed"---        [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Mixed.jacobian (fromPair . bs) [r, s, v, t, k]) 2---        , bench "higherGreeks Double" $ nf (runDouble $ \r s v t k -> (Mixed.hessian (fst . bs) [r, s, v, t, k], Mixed.hessian (snd . bs) [r, s, v, t, k])) 2---        , bench "highererGreeks Double" $ nf (runDouble $ \r s v t k -> Mixed.hessianF (fromPair . bs) [r, s, v, t, k]) 2---        , bench "greeks Float" $ nf (runFloat $ \r s v t k -> Mixed.jacobian (fromPair . bs) [r, s, v, t, k]) 2---        , bench "higherGreeks Float" $ nf (runFloat $ \r s v t k -> (Mixed.hessian (fst . bs) [r, s, v, t, k], Mixed.hessian (snd . bs) [r, s, v, t, k])) 2---        , bench "highererGreeks Float" $ nf (runFloat $ \r s v t k -> Mixed.hessianF (fromPair . bs) [r, s, v, t, k]) 2---        ]+   , bgroup "Mixed"+       [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Mixed.jacobian (fromPair . bs) [r, s, v, t, k]) 2+       , bench "higherGreeks Double" $ nf (runDouble $ \r s v t k -> (Mixed.hessian (fst . bs) [r, s, v, t, k], Mixed.hessian (snd . bs) [r, s, v, t, k])) 2+       , bench "highererGreeks Double" $ nf (runDouble $ \r s v t k -> Mixed.hessianF (fromPair . bs) [r, s, v, t, k]) 2+       , bench "greeks Float" $ nf (runFloat $ \r s v t k -> Mixed.jacobian (fromPair . bs) [r, s, v, t, k]) 2+       , bench "higherGreeks Float" $ nf (runFloat $ \r s v t k -> (Mixed.hessian (fst . bs) [r, s, v, t, k], Mixed.hessian (snd . bs) [r, s, v, t, k])) 2+       , bench "highererGreeks Float" $ nf (runFloat $ \r s v t k -> Mixed.hessianF (fromPair . bs) [r, s, v, t, k]) 2+       ]     ]

@@ -0,0 +1,89 @@+-----------------------------------------------------------------------------+-- |+-- Copyright   :  (c) Edward Kmett 2010-2014+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-- Root finding using Halley's rational method (the second in+-- the class of Householder methods). Assumes the function is three+-- times continuously differentiable and converges cubically when+-- progress can be made.+--+-----------------------------------------------------------------------------++module Numeric.AD.Rank1.Halley+  (+  -- * Halley's Method (Tower AD)+    findZero+  , inverse+  , fixedPoint+  , extremum+  ) where++import Prelude hiding (all)+import Numeric.AD.Internal.Forward (Forward)+import Numeric.AD.Internal.On+import Numeric.AD.Internal.Tower (Tower)+import Numeric.AD.Mode+import Numeric.AD.Rank1.Tower (diffs0)+import Numeric.AD.Rank1.Forward (diff)++-- $setup+-- >>> import Data.Complex++-- | The 'findZero' function finds a zero of a scalar function using+-- Halley'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,1.8571428571428572,1.9997967892704736,1.9999999999994755,2.0]+--+-- >>> last $ take 10 $ findZero ((+1).(^2)) (1 :+ 1)+-- 0.0 :+ 1.0+findZero :: (Fractional a, Eq a) => (Tower a -> Tower a) -> a -> [a]+findZero f = go where+  go x = x : if x == xn then [] else go xn where+    (y:y':y'':_) = diffs0 f x+    xn = x - 2*y*y'/(2*y'*y'-y*y'')+{-# INLINE findZero #-}++-- | The 'inverse' function inverts a scalar function using+-- Halley'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.+--+-- Note: the @take 10 $ inverse sqrt 1 (sqrt 10)@ example that works for Newton's method+-- fails with Halley's method because the preconditions do not hold!+inverse :: (Fractional a, Eq a) => (Tower a -> Tower a) -> a -> a -> [a]+inverse f x0 y = findZero (\x -> f x - auto y) x0+{-# INLINE inverse  #-}++-- | The 'fixedPoint' function find a fixedpoint of a scalar+-- function using Halley'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 :: (Fractional a, Eq a) => (Tower a -> Tower a) -> a -> [a]+fixedPoint f = findZero (\x -> f x - x)+{-# INLINE fixedPoint #-}+++-- | The 'extremum' function finds an extremum of a scalar+-- function using Halley'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.+--+-- >>> take 10 $ extremum cos 1+-- [1.0,0.29616942658570555,4.59979519460002e-3,1.6220740159042513e-8,0.0]+extremum :: (Fractional a, Eq a) => (On (Forward (Tower a)) -> On (Forward (Tower a))) -> a -> [a]+extremum f = findZero (diff (off . f . On))+{-# INLINE extremum #-}