diff --git a/fortran/cobyla/cobylb.f90 b/fortran/cobyla/cobylb.f90 index 8f3ff4d352..484ba32a3c 100644 --- a/fortran/cobyla/cobylb.f90 +++ b/fortran/cobyla/cobylb.f90 @@ -17,7 +17,7 @@ module cobylb_mod ! ! Started: July 2021 ! -! Last Modified: Sunday, March 24, 2024 AM02:53:57 +! Last Modified: Sunday, March 24, 2024 PM05:17:30 !--------------------------------------------------------------------------------------------------! implicit none @@ -40,8 +40,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et ! Common modules use, non_intrinsic :: checkexit_mod, only : checkexit -use, non_intrinsic :: consts_mod, only : RP, IK, ZERO, ONE, HALF, QUART, TENTH, EPS, REALMAX, & - & DEBUGGING, MIN_MAXFILT +use, non_intrinsic :: consts_mod, only : RP, IK, ZERO, ONE, HALF, TENTH, EPS, REALMAX, DEBUGGING, MIN_MAXFILT use, non_intrinsic :: debug_mod, only : assert use, non_intrinsic :: evaluate_mod, only : evaluate use, non_intrinsic :: history_mod, only : savehist, rangehist @@ -132,6 +131,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et real(RP) :: cpen ! Penalty parameter for constraint in merit function (PARMU in Powell's code) real(RP) :: cval(size(x) + 1) real(RP) :: d(size(x)) +real(RP) :: delbar real(RP) :: delta real(RP) :: distsq(size(x) + 1) real(RP) :: dnorm @@ -158,17 +158,6 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et ! 312--314 of Powell's cobylb.f code. Powell's code revises ACTREM to CVAL(N + 1) - CSTRV and PREREM ! to PREREC in this case, which is crucial for feasibility problems. real(RP), parameter :: cpenmin = EPS -! FACTOR_ALPHA, FACTOR_BETA, and FACTOR_GAMMA are factors that COBYLB uses when managing the -! simplex. Note the following. -! 1. FACTOR_ALPHA < FACTOR_GAMMA < 1 < FACTOR_BETA. -! 2. FACTOR_GAMMA has nothing to do with GAMMA1 and GAMMA2, which are the contracting/expanding -! factors for updating the trust-region radius DELTA. -! 3. Powell used one more factor FACTOR_DELTA = 1.1 (in general, 1 < FACTOR_DELTA <= FACTOR_BETA). -! It had nothing to do with DELTA, which is the trust-region radius. It was used when defining -! JDROP_TR. We use a completely different scheme (see SETDROP_TR), which does not need FACTOR_DELTA. -real(RP), parameter :: factor_alpha = QUART ! The factor alpha in the COBYLA paper -real(RP), parameter :: factor_beta = 2.1_RP ! The factor beta in the COBYLA paper -real(RP), parameter :: factor_gamma = HALF ! The factor gamma in the COBYLA paper ! Sizes m_lcon = int(size(bvec), kind(m_lcon)) @@ -202,9 +191,6 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et call assert(size(conhist, 1) == m .and. maxconhist * (maxconhist - maxhist) == 0, & & 'SIZE(CONHIST, 1) == M, SIZE(CONHIST, 2) == 0 or MAXHIST', srname) call assert(maxchist * (maxchist - maxhist) == 0, 'SIZE(CHIST) == 0 or MAXHIST', srname) - call assert(factor_alpha > 0 .and. factor_alpha < factor_gamma .and. factor_gamma < 1, & - & '0 < FACTOR_ALPHA < FACTOR_GAMMA < 1', srname) - call assert(factor_beta > 1, 'FACTOR_BETA > 1', srname) end if !====================! @@ -475,8 +461,8 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et !----------------------------------------------------------------------------------------------! - ! Before the next trust-region iteration, we possibly improve the geometry of simplex or - ! reduces RHO according to IMPROVE_GEO and REDUCE_RHO. Now we decide these indicators. + ! Before the next trust-region iteration, we possibly improve the geometry of simplex or reduce + ! RHO according to IMPROVE_GEO and REDUCE_RHO. Now we decide these indicators. ! N.B.: We must ensure that the algorithm does not set IMPROVE_GEO = TRUE at infinitely many ! consecutive iterations without moving SIM(:, N+1) or reducing RHO. Otherwise, the algorithm ! will get stuck in repetitive invocations of GEOSTEP. This is ensured by the following facts. @@ -570,17 +556,14 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et ! trial point in the current iteration, which was then replaced with the same geometry- ! improving point in the next iteration, and so on. In this process, the simplex alternated ! between two configurations, neither of which had acceptable geometry. Thus RHO was never - ! reduced, leading to infinite cycling. (N.B.: our implementation uses DELTA as the trust + ! reduced, leading to infinite cycling. (N.B.: Our implementation uses DELTA as the trust ! region radius, with RHO being its lower bound. When the infinite cycling occurred in this ! test, DELTA = RHO and it could not be reduced due to the requirement that DELTA >= RHO.) jdrop_geo = int(maxloc(sum(sim(:, 1:n)**2, dim=1), dim=1), kind(jdrop_geo)) ! Calculate the geometry step D. - ! In NEWUOA, GEOSTEP takes DELBAR = MAX(MIN(TENTH * SQRT(MAXVAL(DISTSQ)), HALF * DELTA), RHO) - ! rather than DELTA. This should not be done here, because D should improve the geometry of - ! the simplex when SIM(:, JDROP) is replaced with D; the quality of the geometry is defined - ! by DELTA instead of DELBAR as in (14) of the COBYLA paper. See GEOSTEP for more detail. - d = geostep(jdrop_geo, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamma, simi) + delbar = HALF * delta + d = geostep(jdrop_geo, amat, bvec, conmat, cpen, cval, delbar, fval, simi) x = sim(:, n + 1) + d ! If X is close to one of the points in the interpolation set, then we do not evaluate the @@ -588,8 +571,8 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et ! 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 FACTOR_GAMMA*DELTA, yet X may be close to them due - ! to rounding. In an experiment with single precision on 20240317, X = SIM(:, N+1) occurred. + ! 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. 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 diff --git a/fortran/cobyla/geometry.f90 b/fortran/cobyla/geometry.f90 index a37fdfe8e1..280a66607e 100644 --- a/fortran/cobyla/geometry.f90 +++ b/fortran/cobyla/geometry.f90 @@ -8,7 +8,7 @@ module geometry_cobyla_mod ! ! Started: July 2021 ! -! Last Modified: Saturday, March 23, 2024 PM09:49:28 +! Last Modified: Sunday, March 24, 2024 PM05:13:08 !--------------------------------------------------------------------------------------------------! implicit none @@ -109,7 +109,8 @@ function setdrop_tr(ximproved, d, delta, rho, sim, simi) result(jdrop) !vsig = ONE / sqrt(sum(simi**2, dim=2)) !sigbar = abs(simid) * vsig ! -!! The following JDROP will overwrite the previous one if its premise holds. FACTOR_DELTA = 1.1. +!! The following JDROP will overwrite the previous one if its premise holds. FACTOR_DELTA = 1.1 +!! and FACTOR_ALPHA = 0.25. !mask = (veta > factor_delta * delta .and. (sigbar >= factor_alpha * delta .or. sigbar >= vsig)) !if (any(mask)) then ! jdrop = int(maxloc(veta, mask=mask, dim=1), kind(jdrop)) @@ -194,7 +195,7 @@ function setdrop_tr(ximproved, d, delta, rho, sim, simi) result(jdrop) end function setdrop_tr -function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamma, simi) result(d) +function geostep(jdrop, amat, bvec, conmat, cpen, cval, delbar, fval, simi) result(d) !--------------------------------------------------------------------------------------------------! ! This function calculates a geometry step so that the geometry of the interpolation set is improved ! when SIM(:, JDRO_GEO) is replaced with SIM(:, N+1) + D. See (15)--(17) of the COBYLA paper. @@ -215,8 +216,7 @@ function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamm real(RP), intent(in) :: conmat(:, :) ! CONMAT(M, N+1) real(RP), intent(in) :: cpen real(RP), intent(in) :: cval(:) ! CVAL(N+1) -real(RP), intent(in) :: delta -real(RP), intent(in) :: factor_gamma +real(RP), intent(in) :: delbar real(RP), intent(in) :: fval(:) ! FVAL(N+1) real(RP), intent(in) :: simi(:, :) ! SIMI(N, N) @@ -242,7 +242,7 @@ function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamm if (DEBUGGING) then call assert(m >= m_lcon .and. m >= 0, 'M >= 0', srname) call assert(n >= 1, 'N >= 1', srname) - call assert(delta > 0, 'DELTA > 0', srname) + call assert(delbar > 0, 'DELBAR > 0', srname) call assert(cpen > 0, 'CPEN > 0', srname) call assert(size(simi, 1) == n .and. size(simi, 2) == n, 'SIZE(SIMI) == [N, N]', srname) call assert(all(is_finite(simi)), 'SIMI is finite', srname) @@ -253,7 +253,6 @@ function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamm call assert(size(cval) == n + 1 .and. .not. any(cval < 0 .or. is_nan(cval) .or. is_posinf(cval)), & & 'SIZE(CVAL) == NPT and CVAL does not contain negative NaN/+Inf', srname) call assert(jdrop >= 1 .and. jdrop <= n, '1 <= JDROP <= N', srname) - call assert(factor_gamma > 0 .and. factor_gamma < 1, '0 < FACTOR_GAMMA < 1', srname) end if !====================! @@ -261,12 +260,9 @@ function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamm !====================! ! SIMI(JDROP, :) is a vector perpendicular to the face of the simplex to the opposite of vertex -! JDROP. Set D to the vector in this direction and with length FACTOR_GAMMA * DELTA. Since -! FACTOR_ALPHA < FACTOR_GAMMA < FACTOR_BETA, D improves the geometry of the simplex as per (14) of -! the COBYLA paper. This also explains why this subroutine does not replace DELTA with -! DELBAR = MAX(MIN(TENTH * SQRT(MAXVAL(DISTSQ)), HALF * DELTA), RHO) as in NEWUOA. +! JDROP. Set D to the vector in this direction and with length DELBAR. d = simi(jdrop, :) -d = factor_gamma * delta * (d / norm(d)) +d = delbar * (d / norm(d)) ! The code below chooses the direction of D according to an approximation of the merit function. ! See (17) of the COBYLA paper and line 225 of Powell's cobylb.f. @@ -293,10 +289,10 @@ function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamm ! Postconditions if (DEBUGGING) then call assert(size(d) == n .and. all(is_finite(d)), 'SIZE(D) == N, D is finite', srname) - ! In theory, ||S|| == FACTOR_GAMMA*DELTA, which may be false due to rounding, but not too far. + ! In theory, ||S|| == DELBAR, which may be false due to rounding, but not too far. ! It is crucial to ensure that the geometry step is nonzero, which holds in theory. - call assert(norm(d) > 0.9_RP * factor_gamma * delta .and. norm(d) <= 1.1_RP * factor_gamma * delta, & - & '||D|| == FACTOR_GAMMA*DELTA', srname) + call assert(norm(d) > 0.9_RP * delbar .and. norm(d) <= 1.1_RP * delbar, & + & '||D|| == DELBAR', srname) end if end function geostep