From 084d8a01dcc4b31cb5730ae3b3b0def42b2136e7 Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 18 May 2021 07:47:18 -0400 Subject: [PATCH] Improved exception handling in the BFGS code, and improved convergence by raising the objective function to the 4th power --- src/symba/symba_frag_pos.f90 | 15 +++++++++------ src/util/util_minimize_bfgs.f90 | 27 +++++++++------------------ 2 files changed, 18 insertions(+), 24 deletions(-) diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index e9ed8d3d0..cda9d3f4b 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -543,7 +543,7 @@ subroutine set_fragment_radial_velocities(lmerge) ! Arguments logical, intent(out) :: lmerge ! Internals - real(DP), parameter :: TOL = 1e-6_DP + real(DP), parameter :: TOL = epsilon(1._DP) real(DP), dimension(:), allocatable :: vflat logical :: lerr integer(I4B) :: i @@ -561,7 +561,7 @@ subroutine set_fragment_radial_velocities(lmerge) ! Initialize radial velocity magnitudes with a random value: allocate(v_r_initial, source=v_r_mag) call random_number(v_r_initial(:)) - v_r_initial(:) = v_r_initial(:) * sqrt(2 * ke_target / nfrag / m_frag(:)) + v_r_initial(:) = 3 * v_r_initial(:) * sqrt(2 * ke_target / nfrag / m_frag(:)) ! Initialize the lambda function using a structure constructor that calls the init method ! Minimize error using the BFGS optimizer v_r_mag(:) = util_minimize_bfgs(ke_constraint(ke_objective_function, v_r_unit, v_t_mag, v_t_unit, x_frag, m_frag, L_frag_orb, ke_target), nfrag, v_r_initial, TOL, lerr) @@ -584,7 +584,7 @@ subroutine set_fragment_radial_velocities(lmerge) v_r_mag(i) = dot_product(v_r(:,i), v_r_unit(:, i)) v_frag(:, i) = v_r_mag(i) * v_r_unit(:, i) + v_t_mag(i) * v_t_unit(:, i) end do - !call shift_vector_to_origin(m_frag, v_frag) + call shift_vector_to_origin(m_frag, v_frag) end if do i = 1, nfrag @@ -620,11 +620,14 @@ function ke_objective_function(v_r_mag, v_r_unit, v_t_mag, v_t_unit, x_frag, m_f v_shift(:,i) = v_r_mag(i) * v_r_unit(:, i) end do call shift_vector_to_origin(m_frag, v_shift) - fval = 0.0_DP + + fval = -ke_target do i = 1, nfrag - v_shift(:, i) = v_shift(:, i) + v_t_mag(i) * v_t_unit(:, i) - fval = fval + 0.5 * m_frag(i) * dot_product(v_shift(:, i), v_shift(:, i)) + v_shift(:, i) = v_shift(:, i) + v_t_mag(i) * v_t_unit(:, i) + vcom(:) + fval = fval + 0.5_DP * m_frag(i) * dot_product(v_shift(:, i), v_shift(:, i)) end do + !write(*,*) 'fval: ',fval + fval = fval**2 return diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index fe942109d..ddb715918 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -1,4 +1,4 @@ -function util_minimize_bfgs(f, N, x0, eps1, lerr) result(x1) +function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) !! author: David A. Minton !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! This function implements the Broyden-Fletcher-Goldfarb-Shanno method to determine the minimum of a function of N variables. @@ -21,16 +21,14 @@ function util_minimize_bfgs(f, N, x0, eps1, lerr) result(x1) integer(I4B), intent(in) :: N class(lambda_obj), intent(in) :: f real(DP), dimension(:), intent(in) :: x0 - real(DP), intent(in) :: eps1 + real(DP), intent(in) :: eps logical, intent(out) :: lerr ! Result real(DP), dimension(:), allocatable :: x1 ! Internals integer(I4B) :: i, j, k, l, conv, num, fnum - !integer(I4B), parameter :: MAXLOOP = 10 !! Maximum number of loops before method is determined to have failed - integer(I4B) :: MAXLOOP - !real(DP), parameter :: gradeps = 1e-5_DP !! Tolerance for gradient calculations - real(DP) :: gradeps, eps !! Tolerance for gradient calculations + integer(I4B), parameter :: MAXLOOP = 500 !! Maximum number of loops before method is determined to have failed + real(DP), parameter :: gradeps = 1e-6_DP !! Tolerance for gradient calculations real(DP), dimension(N) :: S !! Direction vectors real(DP), dimension(N) :: Snorm !! normalized direction real(DP), dimension(N,N) :: H !! Approximated inverse Hessian matrix @@ -41,11 +39,6 @@ function util_minimize_bfgs(f, N, x0, eps1, lerr) result(x1) real(DP), dimension(N,N) :: PP, PyH, HyP real(DP) :: yHy, Py - open(unit=42,file='input.txt',status='old') - read(42,*) MAXLOOP - read(42,*) gradeps - read(42,*) eps - close(42) fnum = 0 lerr = .false. ! Initialize approximate Hessian with the identity matrix (i.e. begin with method of steepest descent) @@ -109,14 +102,12 @@ function util_minimize_bfgs(f, N, x0, eps1, lerr) result(x1) end do ! update H matrix H(:,:) = H(:,:) + ((1._DP - yHy / Py) * PP(:,:) - PyH(:,:) - HyP(:,:)) / Py + if (any(H(:,:) > sqrt(huge(1._DP)) / N)) then + write(*,*) 'Did not converge after ',i,'iterations: H too big' + return + end if end do - write(*,*) "Did not converge!" - write(*,*) "Py: ",Py - do k = 1, N - write(*,'("grad(",I0.2,") = ",ES12.4)') k, grad1(k) - end do - write(*,*) "|grad|: Min Mean Max" - write(*,*) minval(abs(grad1(:))), sum(abs(grad1(:))) / N, maxval(abs(grad1(:))) + write(*,*) "Did not converge! Ran out of loops" return contains