Skip to content

Commit

Permalink
Merge pull request #20 from jacobwilliams/17-itp
Browse files Browse the repository at this point in the history 
added ITP method
  • Loading branch information
@jacobwilliams
jacobwilliams authored Sep 4, 2022
2 parents 7a73293 + ee54508 commit f836b63
Show file tree
Hide file tree
Showing 3 changed files with 273 additions and 8 deletions.
3 changes: 2 additions & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@ Procedure | Description
[`brentq`](https://jacobwilliams.github.io/roots-fortran/proc/brentq.html) | SciPy variant of `brent`
[`chandrupatla`](https://jacobwilliams.github.io/roots-fortran/proc/chandrupatla.html) | Hybrid quadratic/bisection algorithm
[`illinois`](https://jacobwilliams.github.io/roots-fortran/proc/illinois.html) | Illinois method
[`itp`](https://jacobwilliams.github.io/roots-fortran/proc/itp.html) | Interpolate Truncate and Project method
[`muller`](https://jacobwilliams.github.io/roots-fortran/proc/muller.html) | Improved Muller method (for real roots only)
[`pegasus`](https://jacobwilliams.github.io/roots-fortran/proc/pegasus.html) | Pegasus method
[`regula_falsi`](https://jacobwilliams.github.io/roots-fortran/proc/regula_falsi.html) | Classic regula falsi method
Expand All @@ -84,7 +85,7 @@ In general, all the methods are guaranteed to converge. Some will be more effici

* Simple classical methods (`bisection`, `regula_falsi`, `illinois`, `ridders`).
* Newfangled methods (`zhang`, `barycentric`, `blendtf`, `bdqrf`, `anderson_bjorck_king`). These rarely or ever seem to be better than the best methods.
* Best methods (`anderson_bjorck`, `muller`, `pegasus`, `toms748`, `brent`, `brentq`, `brenth`, `chandrupatla`). Generally, one of these will be the most efficient method.
* Best methods (`anderson_bjorck`, `muller`, `pegasus`, `toms748`, `brent`, `brentq`, `brenth`, `chandrupatla`, `itp`). Generally, one of these will be the most efficient method.

Note that some of the implementations in this library contain additional checks for robustness, and so may behave better than naive implementations of the same algorithms. In addition, all methods have an option to fall back to bisection if the method fails to converge.

Expand Down
177 changes: 176 additions & 1 deletion src/root_module.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@ module root_module

integer,parameter :: wp = root_module_rk !! local copy of `root_module_rk` with a shorter name
integer,parameter :: name_len = 32 !! max length of the method names
real(wp),parameter,private :: golden_ratio = (1.0_wp + sqrt(5.0_wp))/2.0_wp !! golden ratio `phi`

type,public :: root_method
!! a type to define enums for the different methods
Expand All @@ -67,6 +68,7 @@ module root_module
type(root_method),parameter,public :: root_method_anderson_bjorck_king = root_method(15, 'anderson_bjorck_king') !! enum type for `anderson_bjorck_king` method
type(root_method),parameter,public :: root_method_blendtf = root_method(16, 'blendtf') !! enum type for `blendtf` method
type(root_method),parameter,public :: root_method_barycentric = root_method(17, 'barycentric') !! enum type for `barycentric` method
type(root_method),parameter,public :: root_method_itp = root_method(18, 'itp') !! enum type for `itp` method

type(root_method),parameter,dimension(*),public :: set_of_root_methods = &
[ root_method_brent, &
Expand All @@ -85,7 +87,8 @@ module root_module
root_method_zhang, &
root_method_anderson_bjorck_king, &
root_method_blendtf, &
root_method_barycentric ] !! list of the available methods (see [[root_scalar]])
root_method_barycentric, &
root_method_itp ] !! list of the available methods (see [[root_scalar]])

type,abstract,public :: root_solver
!! abstract class for the root solver methods
Expand Down Expand Up @@ -243,6 +246,21 @@ module root_module
procedure,public :: find_root => barycentric
end type barycentric_solver

type,extends(root_solver),public :: itp_solver
!! ITP root solver
private

! tuning parameters for ITP:
real(wp) :: k1 = 0.1_wp !! from (0, inf)
real(wp) :: k2 = 0.98_wp * (1.0_wp + golden_ratio) !! from [1, 1+phi]
integer :: n0 = 1

contains
private
procedure,public :: find_root => itp
procedure,public :: set_optional_inputs => itp_optional_inputs
end type itp_solver

abstract interface
function func(me,x) result(f)
!! Interface to the function to be minimized
Expand Down Expand Up @@ -437,6 +455,7 @@ subroutine root_scalar_by_type(method,fun,ax,bx,xzero,fzero,iflag,&
case(root_method_anderson_bjorck_king%id); allocate(anderson_bjorck_king_solver :: s)
case(root_method_blendtf%id); allocate(blendtf_solver :: s)
case(root_method_barycentric%id); allocate(barycentric_solver :: s)
case(root_method_itp%id); allocate(itp_solver :: s)

case default
iflag = -999 ! invalid method
Expand Down Expand Up @@ -2450,6 +2469,162 @@ subroutine barycentric(me,ax,bx,fax,fbx,xzero,fzero,iflag)
end subroutine barycentric
!*****************************************************************************************

!*****************************************************************************************
!>
! For the [[itp]] method, set the optional inputs.

subroutine itp_optional_inputs(me,k1,k2,n0)

implicit none

class(itp_solver),intent(inout) :: me

real(wp),intent(in),optional :: k1 !! from (0, inf) [Default is 0.1]
real(wp),intent(in),optional :: k2 !! from [1, 1+phi] [Default is 0.98*(1+phi)]
integer,intent(in),optional :: n0 !! [Default is 1

if (present(k1)) me%k1 = k1
if (present(k2)) me%k2 = k2
if (present(n0)) me%n0 = n0

end subroutine itp_optional_inputs
!*****************************************************************************************

!*****************************************************************************************
!>
! Compute the zero of the function f(x) in the interval ax,bx using the
! Interpolate Truncate and Project (ITP) method.
!
!### See also
! * Oliveira, I. F. D., Takahashi, R. H. C.,
! "An Enhancement of the Bisection Method Average Performance Preserving Minmax Optimality",
! ACM Transactions on Mathematical Software. 47 (1): 5:1-5:24. (2020-12-06)
!
!### Notes
! * This implementation differs from the reference, in that it has an additional
! check to make sure the `a,b` values have changed at each iteration. If they
! haven't, it will perform a basic bisection.

subroutine itp(me,ax,bx,fax,fbx,xzero,fzero,iflag)

implicit none

class(itp_solver),intent(inout) :: me
real(wp),intent(in) :: ax !! left endpoint of initial interval
real(wp),intent(in) :: bx !! right endpoint of initial interval
real(wp),intent(in) :: fax !! `f(ax)`
real(wp),intent(in) :: fbx !! `f(ax)`
real(wp),intent(out) :: xzero !! abscissa approximating a zero of `f` in the interval `ax`,`bx`
real(wp),intent(out) :: fzero !! value of `f` at the root (`f(xzero)`)
integer,intent(out) :: iflag !! status flag (`0`=root found, `-2`=max iterations reached)

real(wp) :: a,b,ya,yb,x12,r,d,xf,bma,sigma,xt,xitp,yitp,term,aprev,bprev,denom
integer :: n12, nmax
integer :: j !! iteration counter
logical :: root_found !! convergence in x
logical :: fail !! if we can't interpolate

real(wp),parameter :: log2 = log(2.0_wp)

! initialize:
if (fax < fbx) then
a = ax
b = bx
ya = fax
yb = fbx
else
a = bx
b = ax
ya = fbx
yb = fax
end if
iflag = 0
term = (b-a)/(2.0_wp*me%rtol)
n12 = ceiling ( log(term) / log2 ) ! ceiling(log2(term))
nmax = n12 + me%n0
aprev = huge(1.0_wp) ! initialize to unusual values
bprev = huge(1.0_wp) !

! main loop
do j = 0, min(nmax,me%maxiter)

x12 = bisect(a,b)
bma = b - a
! note: protect for r<0 as mentioned in paper
r = max(me%rtol * 2.0_wp ** (nmax-j) - bma/2.0_wp, 0.0_wp)
d = me%k1 * bma**me%k2

! interpolation:
denom = yb - ya
fail = abs(denom) <= tiny(1.0_wp) ! check for divide by zero

if (.not. fail) then
xf = (yb*a - ya*b) / denom

! truncation:
sigma = sign(1.0_wp, x12 - xf)
if (d <= abs(x12-xf)) then
xt = xf + sigma * d
else
xt = x12
end if

! projection:
if (abs(xt - x12) <= r) then
xitp = xt
else
xitp = x12 - sigma * r
end if

! updating interval:
yitp = me%f(xitp)
if (solution(xitp,yitp,me%ftol,xzero,fzero)) return
if (yitp > 0.0_wp) then
b = xitp
yb = yitp
elseif (yitp < 0.0_wp) then
a = xitp
ya = yitp
else
a = xitp
b = xitp
end if

end if

if (fail .or. (a==aprev .and. b==bprev)) then
! if the interval hasn't changed, then it is stuck.
! [this can happen in the test cases].
! So just do a bisection
xitp = bisect(a,b)
yitp = me%f(xitp)
if (solution(xitp,yitp,me%ftol,xzero,fzero)) return
if (ya*yitp<0.0_wp) then
! root lies between a and xitp
b = xitp
yb = yitp
else
! root lies between xitp and b
a = xitp
ya = yitp
end if
end if
aprev = a
bprev = b

! check for convergence:
root_found = me%converged(a,b)
if (root_found .or. j==me%maxiter) then
call choose_best(a,b,ya,yb,xzero,fzero)
if (.not. root_found) iflag = -2 ! max iterations reached
exit
end if

end do

end subroutine itp
!*****************************************************************************************

!*****************************************************************************************
!>
! Determines convergence in x based on if the reltol or abstol is satisfied.
Expand Down
101 changes: 95 additions & 6 deletions test/root_tests.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ program root_tests
character(len=*),parameter :: fmt = '( A20,1X,A3,1X,A4,1X,A16, 1X,A25, 1X,A16, 1X,A5,1X,A5,1X,A8 )' !! format for header
character(len=*),parameter :: dfmt = '(1P,A20,1X,I3,1X,I4,1X,E16.6,1X,E25.16,1X,E16.6,1X,I5,1X,I5,1X,E8.1)' !! format for results

integer,parameter :: number_of_methods = 17 !! number of methods to test
integer,parameter :: number_of_methods = 18 !! number of methods to test
character(len=100),dimension(number_of_methods),parameter :: methods = [ &
'bisection ', &
'brent ', &
Expand All @@ -52,7 +52,8 @@ program root_tests
'toms748 ', &
'zhang ', &
'blendtf ', &
'barycentric '] !! method names
'barycentric ', &
'itp '] !! method names

integer,dimension(number_of_methods) :: number_of_wins, ivec, number_of_failures, ivec2

Expand Down Expand Up @@ -707,7 +708,7 @@ subroutine problems(x, ax, bx, fx, xroot, cases, num_of_problems, latex)
root = 0.0_wp
if (present(x)) f = sin(x*10.0_wp) + sin(x*2.0_wp) + &
atan(x/3.0_wp)*x**2 + exp(x) - 1.0_wp
if (present(latex)) latex = '\sin(10 x) + \sin(2 x) + \atan(x/3) x^2 + \mathrm{e}^x - 1'
if (present(latex)) latex = '\sin(10 x) + \sin(2 x) + \tan^{-1} (x/3) x^2 + \mathrm{e}^x - 1'

! functions are from:
! T. R. Chandrupatla, "A new hybrid quadratic/bisection algorithm for finding the zero of
Expand Down Expand Up @@ -1202,21 +1203,109 @@ subroutine problems(x, ax, bx, fx, xroot, cases, num_of_problems, latex)
if (present(x)) f = sin(x) - x**2
if (present(latex)) latex = '\sin(x) - x^2'

case(118) ! https://en.wikipedia.org/wiki/ITP_method#Example:_Finding_the_root_of_a_polynomial
a = 1.0_wp
b = 2.0_wp
root = 1.5213797068045798E+00_wp
if (present(x)) f = x**3 - x - 2.0_wp
if (present(latex)) latex = 'x^3 - x - 2'

! from the ITP paper : Table 1
case(119) ! Lambert
a = -1.0_wp
b = 1.0_wp
root = 5.6714329040978384E-01_wp
if (present(x)) f = x*exp(x) - 1.0_wp
if (present(latex)) latex = 'x \exp(x) - 1'
case(120) ! Trigonometric 1
a = -1.0_wp
b = 1.0_wp
root = 1.0E-01_wp
if (present(x)) f = tan(x - 1.0_wp / 10.0_wp)
if (present(latex)) latex = '\tan( x - \frac{1}{10} )'
case(121) ! Trigonometric 2
a = -1.0_wp
b = 1.0_wp
root = -5.2359877559829893E-01_wp
if (present(x)) f = sin(x) + 1.0_wp / 2.0_wp
if (present(latex)) latex = '\sin(x) + \frac{1}{2}'
case(122) ! Polynomial 1
a = -1.0_wp
b = 1.0_wp
root = -8.4391456864926628E-01_wp
if (present(x)) f = 4.0_wp*x**5 + x**2 + 1.0_wp
if (present(latex)) latex = '4x^5 + x^2 + 1'
case(123) ! Polynomial 2
a = -1.0_wp
b = 1.0_wp
root = 8.3507904272355904E-01_wp
if (present(x)) f = x + x**10 - 1.0_wp
if (present(latex)) latex = 'x + x^{10} - 1'
case(124) ! Exponential
a = -1.0_wp
b = 1.0_wp
root = 8.7356852683023178E-01_wp
if (present(x)) f = pi**x - exp(1.0_wp)
if (present(latex)) latex = '\pi^x - e'
case(125) ! Posynomial
a = -1.0_wp
b = 1.0_wp
root = -3.7020127707860923E-02_wp
if (present(x)) f = 1.0_wp/3.0_wp + sign(1.0_wp, x)*abs(x)**(1.0_wp/3.0_wp) + x**3
if (present(latex)) latex = '1/3 + \sign x |x|^{1/3} + x^3'
case(126) ! Weierstrass
a = -1.0_wp
b = 1.0_wp
root = -2.0065503320580689E-04_wp
if (present(x)) f = 1.0_wp/10.0_wp**3 + sum([(sin(pi*i**3*x/2.0_wp)/pi/i**3, i = 1, 10)])
if (present(latex)) latex = '\frac{1}{10^3} + \sum_{i=1}^{10} \sin (\frac{\pi i^3 x}{2}) / \pi i^3'
case(127) ! Poly.Frac.
a = -1.0_wp
b = 1.0_wp
root = -6.6666666666666666E-01_wp
if (present(x)) f = (x + 2.0_wp/3.0_wp) / (x + 101.0_wp/100.0_wp)
if (present(latex)) latex = '(x + \frac{2}{3}) / (x + \frac{101}{100})'

! ill-behaved functions:
case(128) ! Polynomial 3
a = -1.0_wp
b = 1.0_wp
root = 1.0e-6_wp
if (present(x)) f = (x * 1.0e6_wp - 1.0_wp)**3
if (present(latex)) latex = '(x 10^6 - 1)^3'
case(129) ! Exp.Poly.
a = -1.0_wp
b = 1.0_wp
root = 1.0e-6_wp
if (present(x)) f = exp(x) * (x * 1.0e6_wp - 1.0_wp)**3
if (present(latex)) latex = 'e^x (x 10^6 - 1)^3'
case(130) ! Tan.Poly.
a = -1.0_wp
b = 1.0_wp
root = 3.3333333333333333E-01_wp
if (present(x)) f = (x - 1.0_wp/3.0_wp)**2 * atan( x - 1.0_wp/3.0_wp )
if (present(latex)) latex = '\sign (3x + 1) \left( 1 - \sqrt{1 - (3x+1)^2/9^2} \right)'
case(131) ! Circles
a = -1.0_wp
b = 1.0_wp
root = -3.3333333333333333E-01_wp
if (present(x)) f = sign(1.0_wp,3.0_wp*x+1.0_wp)*(1.0_wp - sqrt(1.0_wp-(3.0_wp*x+1.0_wp)**2/(9**2)))
if (present(latex)) latex = '' !todo

! linear function:
case(118)
case(132)
a = -1000.0_wp
b = 1.0_wp
root = 0.0_wp
if (present(x)) f = x
if (present(latex)) latex = 'x'


case default
write(*,*) 'invalid case: ', nprob
error stop 'invalid case'
end select

if (present(num_of_problems)) num_of_problems = 118
if (present(num_of_problems)) num_of_problems = 132

! outputs:
if (present(ax)) ax = a
Expand Down

0 comments on commit f836b63

Please sign in to comment.