240404.092542.HKT fotran/uobyqa: the same as e7ff2aa, do not evaluate…

… f at a point if it is close to an interpolation point
libprima · Apr 4, 2024 · 411b17b · 411b17b
1 parent 03c90c0
commit 411b17b
Show file tree

Hide file tree

Showing 4 changed files with 66 additions and 39 deletions.
diff --git a/.development b/.development
diff --git a/fortran/cobyla/cobylb.f90 b/fortran/cobyla/cobylb.f90
@@ -17,7 +17,7 @@ module cobylb_mod
 !
 ! Started: July 2021
 !
-! Last Modified: Sunday, March 31, 2024 PM07:43:03
+! Last Modified: Wednesday, April 03, 2024 PM08:37:51
 !--------------------------------------------------------------------------------------------------!
 
 implicit none
@@ -347,7 +347,6 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
     ! Calculate the trust-region trial step D. Note that D does NOT depend on CPEN.
     d = trstlp(A, -conmat(:, n + 1), delta, g)
     dnorm = min(delta, norm(d))
-    x = sim(:, n + 1) + d
 
     ! Is the trust-region trial step short? Note that we compare DNORM with RHO, not DELTA.
     ! Powell's code essentially defines SHORTD by SHORTD = (DNORM < HALF * RHO). In our tests,
@@ -379,9 +378,11 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
             delta = rho  ! Set DELTA to RHO when it is close to or below.
         end if
     else
+        ! Calculate the next value of the objective and constraint functions.
         ! If X is close to one of the points in the interpolation set, then we do not evaluate the
         ! objective and constraints at X, assuming them to have the values at the closest point.
         ! N.B.: If this happens, do NOT include X into the filter, as F and CONSTR are inaccurate.
+        x = sim(:, n + 1) + d
         distsq(n + 1) = sum((x - sim(:, n + 1))**2)
         distsq(1:n) = [(sum((x - (sim(:, n + 1) + sim(:, j)))**2), j=1, n)]  ! Implied do-loop
         !!MATLAB: distsq(1:n) = sum((x - (sim(:,1:n) + sim(:, n+1)))**2, 1)  % Implicit expansion
@@ -564,15 +565,16 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
         ! Calculate the geometry step D.
         delbar = HALF * delta
         d = geostep(jdrop_geo, amat, bvec, conmat, cpen, cval, delbar, fval, simi)
-        x = sim(:, n + 1) + d
 
+        ! Calculate the next value of the objective and constraint functions.
         ! If X is close to one of the points in the interpolation set, then we do not evaluate the
         ! objective and constraints at X, assuming them to have the values at the closest point.
         ! N.B.:
         ! 1. If this happens, do NOT include X into the filter, as F and CONSTR are inaccurate.
         ! 2. In precise arithmetic, the geometry improving step ensures that the distance between X
         ! and any interpolation point is at least DELBAR, yet X may be close to them due to
         ! rounding. In an experiment with single precision on 20240317, X = SIM(:, N+1) occurred.
+        x = sim(:, n + 1) + d
         distsq(n + 1) = sum((x - sim(:, n + 1))**2)
         distsq(1:n) = [(sum((x - (sim(:, n + 1) + sim(:, j)))**2), j=1, n)]  ! Implied do-loop
         !!MATLAB: distsq(1:n) = sum((x - (sim(:,1:n) + sim(:, n+1)))**2, 1)  % Implicit expansion
@@ -645,7 +647,8 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
 end do  ! End of DO TR = 1, MAXTR. The iterative procedure ends.
 
 ! Return from the calculation, after trying the last trust-region step if it has not been tried yet.
-! Ensure that X has not been updated after SHORTD == TRUE occurred, or the code below is incorrect.
+! Ensure that D has not been updated after SHORTD == TRUE occurred, or the code below is incorrect.
+x = sim(:, n + 1) + d
 if (info == SMALL_TR_RADIUS .and. shortd .and. norm(x - sim(:, n + 1)) > 1.0E-3_RP * rhoend .and. nf < maxfun) then
     call evaluate(calcfc_internal, x, f, constr)
     cstrv = maximum([ZERO, constr])

diff --git a/fortran/common/history.f90 b/fortran/common/history.f90
@@ -7,7 +7,7 @@ module history_mod
 !
 ! Started: July 2020
 !
-! Last Modified: Sunday, March 24, 2024 AM03:01:46
+! Last Modified: Thursday, April 04, 2024 AM09:21:08
 !--------------------------------------------------------------------------------------------------!
 
 implicit none

diff --git a/fortran/uobyqa/uobyqb.f90 b/fortran/uobyqa/uobyqb.f90
@@ -8,7 +8,7 @@ module uobyqb_mod
 !
 ! Started: February 2022
 !
-! Last Modified: Wednesday, April 03, 2024 AM12:45:05
+! Last Modified: Thursday, April 04, 2024 AM01:10:44
 !--------------------------------------------------------------------------------------------------!
 
 implicit none
@@ -94,6 +94,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r
 ! Local variables
 character(len=*), parameter :: solver = 'UOBYQA'
 character(len=*), parameter :: srname = 'UOBYQB'
+integer(IK) :: k
 integer(IK) :: knew_geo
 integer(IK) :: knew_tr
 integer(IK) :: kopt
@@ -284,9 +285,9 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r
     g = pq(1:n) + smat_mul_vec(pq(n + 1:npt - 1), xpt(:, kopt))
     h = vec2smat(pq(n + 1:npt - 1))
     call trstep(delta, g, h, trtol, d, crvmin)
+    dnorm = min(delta, norm(d))
 
     ! Check whether D is too short to invoke a function evaluation.
-    dnorm = min(delta, norm(d))
     shortd = (dnorm <= HALF * rho)  ! `<=` works better than `<` in case of underflow.
 
     ! Set QRED to the reduction of the quadratic model when the move D is made from XOPT. QRED
@@ -302,21 +303,24 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r
         end if
     else
         ! Calculate the next value of the objective function.
+        ! If X is close to one of the points in the interpolation set, then we do not evaluate the
+        ! objective function X, assuming it to have the value at the closest point.
         x = xbase + (xpt(:, kopt) + d)
-        call evaluate(calfun, x, f)
-        nf = nf + 1_IK
+        distsq = [(sum((x - (xbase + xpt(:, k)))**2, dim=1), k=1, npt)]  ! Implied do-loop
+        !!MATLAB: distsq = sum((x - (xbase + xpt))**2, 1)  % Implicit expansion
+        if (any(distsq <= (1.0E-4 * rhoend)**2)) then
+            k = int(minloc(distsq, mask=(.not. is_nan(distsq)), dim=1), kind(k))
+            f = fval(k)
+        else
+            ! Evaluate the objective function at X, taking care of possible Inf/NaN values.
+            call evaluate(calfun, x, f)
+            nf = nf + 1_IK
+            ! Save X and F into the history.
+            call savehist(nf, x, xhist, f, fhist)
+        end if
 
         ! Print a message about the function evaluation according to IPRINT.
         call fmsg(solver, 'Trust region', iprint, nf, delta, f, x)
-        ! Save X, F into the history.
-        call savehist(nf, x, xhist, f, fhist)
-
-        ! Check whether to exit
-        subinfo = checkexit(maxfun, nf, f, ftarget, x)
-        if (subinfo /= INFO_DFT) then
-            info = subinfo
-            exit
-        end if
 
         ! Update DNORM_REC and MODERR_REC.
         ! DNORM_REC records the DNORM of the recent function evaluations with the current RHO.
@@ -353,6 +357,13 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r
             end if
             ddmove = sum((xdrop - xpt(:, kopt))**2)  ! KOPT is updated.
         end if
+
+        ! Check whether to exit
+        subinfo = checkexit(maxfun, nf, f, ftarget, x)
+        if (subinfo /= INFO_DFT) then
+            info = subinfo
+            exit
+        end if
     end if  ! End of IF (SHORTD .OR. TRFAIL). The normal trust-region calculation ends.
 
 
@@ -456,21 +467,24 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r
         d = geostep(knew_geo, kopt, delbar, pl, xpt)
 
         ! Calculate the next value of the objective function.
+        ! If X is close to one of the points in the interpolation set, then we do not evaluate the
+        ! objective function X, assuming it to have the value at the closest point.
         x = xbase + (xpt(:, kopt) + d)
-        call evaluate(calfun, x, f)
-        nf = nf + 1_IK
+        distsq = [(sum((x - (xbase + xpt(:, k)))**2, dim=1), k=1, npt)]  ! Implied do-loop
+        !!MATLAB: distsq = sum((x - (xbase + xpt))**2, 1)  % Implicit expansion
+        if (any(distsq <= (1.0E-4 * rhoend)**2)) then
+            k = int(minloc(distsq, mask=(.not. is_nan(distsq)), dim=1), kind(k))
+            f = fval(k)
+        else
+            ! Evaluate the objective function at X, taking care of possible Inf/NaN values.
+            call evaluate(calfun, x, f)
+            nf = nf + 1_IK
+            ! Save X and F into the history.
+            call savehist(nf, x, xhist, f, fhist)
+        end if
 
         ! Print a message about the function evaluation according to IPRINT.
         call fmsg(solver, 'Geometry', iprint, nf, delbar, f, x)
-        ! Save X, F into the history.
-        call savehist(nf, x, xhist, f, fhist)
-
-        ! Check whether to exit
-        subinfo = checkexit(maxfun, nf, f, ftarget, x)
-        if (subinfo /= INFO_DFT) then
-            info = subinfo
-            exit
-        end if
 
         ! Update DNORM_REC and MODERR_REC.
         ! DNORM_REC records the DNORM of the recent function evaluations with the current RHO.
@@ -487,6 +501,13 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r
             info = NAN_INF_MODEL
             exit
         end if
+
+        ! Check whether to exit
+        subinfo = checkexit(maxfun, nf, f, ftarget, x)
+        if (subinfo /= INFO_DFT) then
+            info = subinfo
+            exit
+        end if
     end if  ! End of IF (IMPROVE_GEO). The procedure of improving geometry ends.
 
     ! The calculations with the current RHO are complete. Enhance the resolution of the algorithm
@@ -528,22 +549,25 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r
 deallocate (pl)
 
 ! Return from the calculation, after trying the Newton-Raphson step if it has not been tried yet.
-if (info == SMALL_TR_RADIUS .and. shortd .and. dnorm >= TENTH * rhoend .and. nf < maxfun) then
-    x = xbase + (xpt(:, kopt) + d)
+! Ensure that D has not been updated after SHORTD == TRUE occurred, or the code below is incorrect.
+x = xbase + (xpt(:, kopt) + d)
+if (info == SMALL_TR_RADIUS .and. shortd .and. norm(x - (xbase + xpt(:, kopt))) >= TENTH * rhoend .and. nf < maxfun) then
     call evaluate(calfun, x, f)
     nf = nf + 1_IK
+    ! Save X, F into the history.
+    call savehist(nf, x, xhist, f, fhist)
     ! Print a message about the function evaluation according to IPRINT.
     ! Zaikun 20230512: DELTA has been updated. RHO is only indicative here. TO BE IMPROVED.
     call fmsg(solver, 'Trust region', iprint, nf, rho, f, x)
-    ! Save X, F into the history.
-    call savehist(nf, x, xhist, f, fhist)
+    if (f < fval(kopt)) then
+        xpt(:, kopt) = xpt(:, kopt) + d
+        fval(kopt) = f
+    end if
 end if
 
-! Choose the [X, F] to return: either the current [X, F] or [XBASE + XOPT, FOPT].
-if (fval(kopt) < f .or. is_nan(f)) then
-    x = xbase + xpt(:, kopt)
-    f = fval(kopt)
-end if
+! Choose the [X, F] to return.
+x = xbase + xpt(:, kopt)
+f = fval(kopt)
 
 ! Arrange FHIST and XHIST so that they are in the chronological order.
 call rangehist(nf, xhist, fhist)