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

@@ -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
 

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