program pnelmead ! test problem for nelmead -- general minimizer of functions ! IMPLICIT NONE REAL(KIND=2), EXTERNAL :: f0,f1,f2,f3,f4 REAL(KIND=2), DIMENSION(10) :: x,typx REAL(KIND=2) fnew integer i,n,termcd ! data x/ -1.2d0, 1.d0, -1.2d0, 1.d0, 6*0.d0 / data typx/ 10*1.d0 / ! 21 format(i4,f12.6) 22 format(6f13.6) 24 format(//'Zero Test of nelmead -- quadratic function, n=2') 25 format(//'First Test of nelmead -- Rosenbrock function, n=2') 26 format(//'Second Test of nelmead -- Rosenbrock function, n=4') 27 format(//'Third Test of nelmead -- Cox likelihood, n=2') 28 format(//'Fourth Test of nelmead -- D&S Example 9.4.1') ! OPEN( UNIT=6, FILE='nelmead.out' ) ! zero problem in 2 dimensions write(6,24) n = 2 call nelmead(f0,x,n,80,typx,0.1d0,1.d-8,fnew,6,termcd) write(6,21) termcd,fnew write(6,22) (x(i),i=1,n) ! first problem in 2 dimensions write(6,25) n = 2 x(1) = -1.2 x(2) = 1.0 call nelmead(f1,x,n,100,typx,0.1d0,1.d-8,fnew,6,termcd) write(6,21) termcd,fnew write(6,22) (x(i),i=1,n) ! second problem in 4 dimensions write(6,26) n = 4 ! x(1) = -1.2 x(2) = 1.0 x(3) = -1.2 x(4) = 1.0 call nelmead(f2,x,n,750,typx,0.1d0,1.d-8,fnew,6,termcd) write(6,21) termcd,fnew write(6,22) (x(i),i=1,n) ! third problem is Cox loglikelihood write(6,27) n = 2 x(1) = 0.d0 x(2) = 0.d0 call nelmead(f3,x,n,80,typx,0.1d0,1.d-8,fnew,6,termcd) write(6,21) termcd,fnew write(6,22) (x(i),i=1,n) ! fourth problem write(6,28) n = 2 x(1) = 0.d0 x(2) = 0.d0 call nelmead(f4,x,n,80,typx,0.1d0,1.d-8,fnew,6,termcd) write(6,21) termcd,fnew write(6,22) (x(i),i=1,n) stop end program pnelmead REAL(KIND=2) function f0(x) implicit none REAL(KIND=2), DIMENSION(4), INTENT(IN) :: x ! quadratic function in 2 dimensions f0 = 8.*(x(1)-1.)*(x(1)-1.) -4.*(x(1)-1.)*(x(2)-2.) & & + 5.*(x(2)-2.)*(x(2)-2.) + 3. write(6,21) f0,x(1),x(2) 21 format(3f12.6) return end function f0 REAL(KIND=2) function f1(x) implicit none REAL(KIND=2), DIMENSION(4), INTENT(IN) :: x ! rosenbrock function in 2 dimensions f1 = 100.d0*(x(2)-x(1)*x(1))*(x(2)-x(1)*x(1)) & & + (1.d0-x(1))*(1.d0-x(1)) write(6,21) f1,x(1),x(2) 21 format(f16.4,4f12.6) return end function f1 REAL(KIND=2) function f2(x) implicit none REAL(KIND=2), DIMENSION(10), INTENT(IN) :: x integer i,i2,i2m ! extended rosenbrock function in 4 dimensions f2 = 0.d0 do i = 1,2 i2 = i + i i2m = i2 - 1 f2 = f2 + 100.d0*(x(i2)-x(i2m)*x(i2m))*(x(i2)-x(i2m)*x(i2m)) & & + (1.d0-x(i2m))*(1.d0-x(i2m)) end do ! loop on i return end function f2 REAL(KIND=2) function f3(x) IMPLICIT none REAL(KIND=2), DIMENSION(2), INTENT(IN) :: x REAL(KIND=2), DIMENSION(4) :: t REAL(KIND=2) v INTEGER, DIMENSION(4) :: n,s integer i ! Cox logistic regression example data n/ 55, 157, 159, 16 / data s/ 0, 2, 7, 3 /, t/ 7.d0, 14.d0, 27.d0, 51.d0 / f3 = 0.d0 do i = 1,4 v = x(1) + x(2)*t(i) f3 = f3 - v*s(i) + n(i)*dlog(1.d0+dexp(v) ) end do write(6,21) f3,x(1),x(2) 21 format(f16.4,4f12.6) return end function f3 REAL(KIND=2) function f4(x) implicit none REAL(KIND=2), DIMENSION(2), INTENT(IN) :: x ! D&S Example 9.4.1 f4 = (x(1)-2.)**4 + (x(1)-2.d0)*(x(1)-2.d0)*x(2)*x(2) & & + (x(2)+1.d0)*(x(2)+1.d0) write(6,21) f4,x(1),x(2) 21 format(f16.4,4f12.6) return end function f4 subroutine nelmead(f,x,n,mit,typx,scale,epsf,fnew,outcod,termcd) ! subroutine to minimize f ! requiring only function evaluations f(x) ! following simplex algorithm of JA Nelder and R Mead ! The Computer Journal, Volume 7 (1965) pp. 308-313 ! ! f function to be minimized (function subprogram) ! x real vector in: starting point for search ! out: minimum ! n integer dimension of x ! mit integer maximum number of iterations ! typx real vector typical values for |x|, used for scaling ! scale real algorithm scale factor sets relative step size ! 1/8 typical, smaller means smaller steps ! epsf real stop criterion on changes of f in simplex ! 1.d-6 typical, smaller means more work ! fnew real out: value of function at minimum ! outcod integer output code for writing messages ! termcd integer termination code ! termcd=1 means stop -- f changes too small ! termcd=2 means stop -- too many iterations ! ! j f monahan (june, 2005) ! implicit none REAL(KIND=2) f INTEGER, INTENT(IN) :: N,MIT REAL(KIND=2), DIMENSION(N), INTENT(IN OUT) :: x,typx real(kind=2), intent(in) :: scale, epsf REAL(KIND=2), INTENT(IN OUT) :: fnew INTEGER, INTENT(IN) :: outcod INTEGER, INTENT(OUT) :: termcd REAL(KIND=2), DIMENSION( N ) :: xnew, xnewer, xbar, xcur REAL(KIND=2), DIMENSION( N+1 ) :: ylist REAL(KIND=2), DIMENSION( N, N+1 ) :: simplex REAL(KIND=2) ylow,yhigh,ynew,ynewer,y2nd,vary ! integer i,j,kit,high,low,np1 ! ! 21 format(' begin iteration ',i4,' out of max ',i4,' current f'& & ,f10.4,' current x '/2x,6f12.6/4x,6f12.6/4x,6f12.6/6x,6f12.6) 22 format(' expand: ynewer, ynew, ylow ',3f10.4) 23 format(' roll xnew ') 24 format(' xnew worst remaining -- contract') 31 format(' variation of y',f12.6) 41 format(' stopped -- function values not changing over simplex') 42 format(' stopped -- too many iterations') 51 format(15x,i5,6f12.6/6f12.6) ! ! initialization of simplex, etc np1 = n + 1 ! simplest form of simplex simplex = 0. do i = 1, n simplex(i,i) = 1. end do ! loop on i ! now relocate, rescale, and evaluate xbar = 0. do j = 1, np1 do i = 1, n xcur(i) = x(i) + scale*typx(i)*simplex(i,j) simplex(i,j) = xcur(i) xbar(i) = xbar(i) + simplex(i,j) end do ! loop on i ylist(j) = f(xcur) end do ! loop on j xbar = xbar/real(n) !(1) has all pts, (2) divide by n low = np1 ! just for first write ! ! start of main loop ! do kit = 1, mit write(outcod,21) kit,mit,ylow,(simplex(i,low),i=1,n) ! ! find high, low, and should we stop? high = 1 yhigh = ylist(1) low = 1 ylow = ylist(1) do j = 2, np1 if( ylist(j) .lt. ylow ) then low = j ylow = ylist(j) endif ! ( ylist(j) .lt. ylow ) if( ylist(j) .gt. yhigh ) then high = j yhigh = ylist(j) endif ! ( ylist(j) .gt. yhigh ) end do ! loop on j ! measure variation in points vary = (yhigh-ylow) / (1. + abs(ylow)) ! relative error write(outcod, 31) vary if( vary .le. epsf ) then ! stop if yvalues too close ! copy solution do i = 1, n x(i) = simplex(i,low) end do ! loop on i fnew = ylow write(outcod,41) termcd = 1 return ! we're done end if ! ( vary .le. epsf ) ! start to work ! ! first fix xbar by deleting xhigh (that's why we divide by n) do i = 1, n xbar(i) = xbar(i) - simplex(i,high)/real(n) end do ! loop on i ! reflect from xhigh do i = 1,n xnew(i) = 2.*xbar(i) - simplex(i,high) end do ! loop on i ! evaluate ynew = f(xnew) ! first big branch if( ynew .lt. ylow ) then ! expand? do i = 1,n xnewer(i) = 3.*xnew(i) - 2.*xbar(i) end do ! loop on i ! evaluate ynewer = f(xnewer) ! nested branch if ( ynewer .lt. ylow ) then write(outcod,22) ynewer,ynew,ylow ! replace high with xnewer & go do i = 1, n xbar(i) = xbar(i) + xnewer(i)/real(n) simplex(i,high) = xnewer(i) end do ! loop on i ylist(high) = ynewer else ! ( ynewer .lt. ylow ) write(outcod,23) ! replace high with xnew & go do i = 1, n xbar(i) = xbar(i) + xnew(i)/real(n) simplex(i,high) = xnew(i) end do ! loop on i ylist(high) = ynew end if ! ( ynewer .lt. ylow ) ! it's not the best, it is not the worst? else ! ( ynew .lt. ylow ) ! compare to others -- is it better than ! anybody but yhigh? y2nd = ylow do j = 1, np1 if( j .ne. high ) y2nd = max( y2nd, ylist(j) ) end do ! loop on j ! y2nd now second largest if( ynew .le. y2nd ) then write(outcod,23) ! replace high with xnew & go do i = 1, n xbar(i) = xbar(i) + xnew(i)/real(n) simplex(i,high) = xnew(i) end do ! loop on i ylist(high) = ynew else ! ( ynew .le. y2nd ) write(outcod,24) ! complicated branches -- new point is worst! if( ynew .le. yhigh ) then do i = 1, n simplex(i,high) = xnew(i) end do ! loop on i ylist(high) = ynew end if ! ( ynew .le. yhigh ) ! get newer point do i = 1, n xnewer(i) = (xbar(i)+simplex(i,high))/2. end do ! loop on i ! evaluate ynewer = f(xnewer) ! now what to do? if( ynewer .gt. ylist(high) ) then ! replace all points with contraction ! and evaluate xbar = 0. do j = 1, np1 do i = 1, n simplex(i,j) = (simplex(i,j)+simplex(i,low))/2. xbar(i) = xbar(i) + simplex(i,j) end do ! loop on i end do ! loop on j xbar = xbar / real(n) ! evaluate new simplex low = 1 do j = 1, np1 do i = 1, n xcur(i) = simplex(i,j) end do ! loop on i ylist(j) = f(xcur) if( ylist(j) .le. ylist(low) ) low = j end do ! loop on j ylow = ylist(low) ! if we're maxit -- finish with good pt else ! ( ynewer .gt. ylist(high) ) do i = 1, n xbar(i) = xbar(i) + xnewer(i)/real(n) simplex(i,high) = xnewer(i) end do ! loop on i ylist(high) = ynewer end if ! ( ynewer .gt. ylist(high) ) end if ! ( ynew .le. y2nd ) end if ! ( ynew .lt. ylow ) end do ! main loop on kit ! stop because of too many iterations write(outcod,42) termcd = 2 ! copy solution do i = 1, n x(i) = simplex(i,low) end do ! loop on i fnew = ylow return end subroutine nelmead