{-| 

    Module      :  Data.Number.ER.Real.Arithmetic.Taylor
    Description :  interval Newton method
    Copyright   :  (c) Amin Farjudian, Michal Konecny
    License     :  BSD3

    Maintainer  :  mik@konecny.aow.cz
    Stability   :  alpha
    Portability :  portable

    Interval Newton's method for root finding. 
       
    To be used for obtaining functions out of their inverse(s) over various 
    intervals.
-}
module Data.Number.ER.Real.Arithmetic.Newton 
(
    erNewton_FullArgs,
    erNewton_mdfd_FullArgs
)
where

import qualified Data.Number.ER.Real.Approx as RA
import Data.Number.ER.BasicTypes
import Data.Number.ER.Real.Arithmetic.Taylor

erNewton_FullArgs
	:: (RA.ERIntApprox ira)
	=> (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) -- ^ a function and its derivative
	-> ira -- ^ a starting point
	-> ira -- ^ a lower bound of the absolute value of the derivative over the working interval
	-> Int -- ^ number of iterations
	-> EffortIndex  -- ^ the initial index to use for argument function and its derivative
	-> ira -- ^ the result
	
erNewton_FullArgs (f ,df) startPt minDrv iterCnt i = 
    erNewton_FullArgs_aux startPt startOtherPt iterCnt
	where 
    erNewton_FullArgs_aux newtonPt otherPt iterCnt
        | (iterCnt <= 0 || RA.getPrecision result >= prec) =
            result 
        | otherwise = 
            erNewton_FullArgs_aux newNewtonPt newOtherPt (iterCnt - 1)
        where 
        result = 
            newtonPt RA.\/ otherPt
        prec = effIx2prec i 
        newNewtonPt = 
            midPoint $ RA.bounds $ 
            (newtonPt - ( (f i newtonPt) / (( df i newtonPt)))) 
               --  /\ (ira2ra ((ra2ira minDrv) \/ 100000000)))))
        newOtherPt = otherEndPoint newNewtonPt
    startOtherPt = otherEndPoint startPt
    otherEndPoint a =  a - ((f i a) / minDrv) --   /\ (0 \/ 10000000)

    
{-|
    This auxiliary function returns the average of two ra's
-}
midPoint
    :: (RA.ERIntApprox ira)
    => (ira ,ira)
    -> ira
midPoint (x1, x2) = (x1 + x2) / 2
        

{-| Modified Newton Method
    Notes:
    
        1. It has a cubic convergence speed, as opposed to the original Newton's
            square convergence speed.
            
        2. It does not deal with multiple roots.
        
        3. Per iteration, it makes two queries on the derivative, so it best 
            suits the cases where computation of the derivative is at most as
            expensive as the function itself.
-}
erNewton_mdfd_FullArgs
    :: (RA.ERIntApprox ira)
    => (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) -- ^ a function and its derivative
    -> ira -- ^ a starting point
    -> ira -- ^ The minimum of absolute value of derivative over the working interval
    -> Int -- ^ number of iterations
    -> EffortIndex  -- ^ It triggers the initial index to be called by the argument function and its derivative.
    -> ira -- ^ the result
    
erNewton_mdfd_FullArgs (f ,df) startPt minDrv iterCnt i = 
    erNewton_FullArgs_aux startPt startOtherPt iterCnt
    where 
    erNewton_FullArgs_aux newtonPt otherPt iterCnt
        | iterCnt <= 0 = newtonPt RA.\/ otherPt
        | otherwise = erNewton_FullArgs_aux newNewtonPt newOtherPt (iterCnt - 1)
        where
        valueAtNewOtherPt = f i newOtherPt
        derivAtNewtonPt   = df i newOtherPt
        unblurredDeriv = midPoint $ RA.bounds $ derivAtNewtonPt
        intermediaryPt = midPoint $ RA.bounds $ newtonPt - valueAtNewOtherPt / (2 * derivAtNewtonPt)
        derivAtIntermediaryPt = df i intermediaryPt
        newNewtonPt = 
            midPoint $ RA.bounds $ 
            (newtonPt - ( valueAtNewOtherPt / derivAtIntermediaryPt))
        newOtherPt = otherEndPoint newNewtonPt
    startOtherPt = otherEndPoint startPt
    otherEndPoint a = a - ((f i a) / minDrv)

erNewton_mdfd
    :: (RA.ERIntApprox ira)
    => (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) -- ^ a function and its derivative
    -> ira -- ^ a starting point
    -> ira -- ^ The minimum of absolute value of derivative over the working interval
    -> EffortIndex  -- ^ It triggers the initial index to be called by the argument function and its derivative.
    -> ira -- ^ the result
    
erNewton_mdfd (f ,df) startPt minDrv i = 
    erNewton_mdfd_FullArgs (f, df) startPt minDrv (fromInteger $ toInteger $ i) i

--apNewton_mdfd
--    :: (RA.ERIntApprox ira)
--    => (EffortIndex -> ra -> ra, EffortIndex -> ra -> ra) -- ^ a function and its derivative
--    -> ra -- ^ a starting point
--    -> ra -- ^ The minimum of absolute value of derivative over the working interval
--    -> EffortIndex  -- ^ It triggers the initial index to be called by the argument function and its derivative. Moreover, the number of iterations are predefined by this argument.
--    -> ra -- ^ the result
--    
--apNewton_mdfd (f, df) startPt minDrv i =
--    erNewton_mdfd_FullArgs
--
			
--id_RA 
--	:: (RA.ERIntApprox ira)
--	=> EffortIndex -> ira -> ira
--
--id_RA i x = x
--
--const_one_RA
--	:: (RA.ERIntApprox ira)
--	=> EffortIndex -> ira -> ira
--
--const_one_RA i x = (setMinGranularity (effIx2gran i) 1)
--	
--
--test_erNewton_FullArgs_01_RA 
--	:: (RA.ERIntApprox ira)
--	=> EffortIndex -> ira -> ira
--
--test_erNewton_FullArgs_01_RA i x = erNewton_FullArgs_01 ( id_RA, const_one_RA) x 10 i
--
--test_erNewton_FullArgs_01
--	:: (RA.ERIntApprox ira)
--	=> (ConvergRealSeq ira) -> (ConvergRealSeq ira)
--	
--test_erNewton_FullArgs_01 = convertFuncRA2Seq test_erNewton_FullArgs_01_RA
--
--exp_Ra_minus_r_RA
--	:: (RA.ERIntApprox ira)
--	=> EffortIndex -> ira -> ira -> ira 
--	
--exp_Ra_minus_r_RA i r x = (erExp_RA i x) - r
--
--exp_Ra_minus_r 
--	:: (RA.ERIntApprox ira)
--	=> (ConvergRealSeq ira) -> (ConvergRealSeq ira) -> (ConvergRealSeq ira)
--
--exp_Ra_minus_r = convertBinFuncRA2Seq exp_Ra_minus_r_RA
--
--logNewton_RA_02
--    :: (RA.ERIntApprox ira)
--    => EffortIndex -> ira -> ira
--    
--logNewton_RA_02 i x = 
--    erNewton_FullArgs_02
--        ( \ i y -> (erExp_RA i y) - x, erExp_RA) 
--        (setMinGranularity (effIx2gran i) (2)) 
--        (setMinGranularity (effIx2gran i) 1) 
--        i   
--
--logNewton_02  
--    :: (RA.ERIntApprox ira)
--    => (ConvergRealSeq ira) -> (ConvergRealSeq ira)
--    
--logNewton_02 = convertFuncRA2Seq logNewton_RA_02


--logNewton_mdfd_RA
--    :: (RA.ERIntApprox ira)
--    => EffortIndex -> ira -> ira
    
--logNewton_mdfd_RA i x = 
--    erNewton_mdfd_FullArgs
--        ( \ i y -> (erExp_RA i y) - x, erExp_RA) 
--        (setMinGranularity (effIx2gran i) (2)) 
--        (setMinGranularity (effIx2gran i) 1) 
--        i   

--logNewton_mdfd
--    :: (RA.ERIntApprox ira)
--    => (ConvergRealSeq ira) -> (ConvergRealSeq ira)
--    
--logNewton_mdfd = convertFuncRA2Seq logNewton_mdfd_RA