minimizer_tests.F90

Go to the documentation of this file.

!! Copyright (C) 2024 N. Tancogne-Dejean
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module minimizer_tests_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use iso_c_binding
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use parser_oct_m
  use profiling_oct_m
  use messages_oct_m
  use minimizer_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                    &
    test_optimizers
 
  character(len=:), allocatable :: current_log
 
contains
 
  !---------------------------------------------------------------------------
  subroutine test_optimizers(namespace)
    type(namespace_t), intent(in) :: namespace
 
    real(real64), dimension(2) :: x0, x, mass
    real(real64) :: tol_grad = 1.0e-6_real64
    real(real64) :: dt, energy
    integer :: ierr, integrator
    integer, parameter :: max_iter = 1000
 
    push_sub(test_optimizers)
 
    ! See src/main/geom_opt.F90 for the variable description
    call parse_variable(namespace, 'GOFireIntegrator', option__gofireintegrator__verlet, integrator)
 
    !=== Rosenbrock =========================================================
    call set_log_file("rosenbrock.log")
    call io_rm(trim(current_log), namespace)
 
    ! Testing the FIRE algorithm with the Rosenbrock function in 2D
    ! Starting point
    x0 = [3.0_real64, 4.0_real64]
 
    ! Initial parameters for the FIRE algorithm with Velocity Verlet
    x = x0
    ! This timestep is taken to give a reasonable result. Too large vaues lead to unstable results
    ! This is set such that run this through valgrind+verrou gives reproducible results
    dt = 0.001_real64
    mass = m_one
    call minimize_fire(2, 2, x, dt, tol_grad, max_iter, rosenbrock_gradient_2d, write_iter_info, energy, ierr, mass, integrator)
 
    write(message(1), "(a,i6,a)") "FIRE (Rosenbrock) algorithm converged in ", -ierr, " iterations"
    write(message(2), "(a,2(f8.4,a))") " The minimum is found to be at (", x(1), ", ", x(2), ")"
    write(message(3), "(a)") " Analytical minimum is (1,1)."
    call messages_info(3)
    call messages_new_line()
 
    !=== Sphere =============================================================
    call set_log_file("sphere.log")
    call io_rm(trim(current_log), namespace)
 
    x0 = [3.0_real64, -4.0_real64]
    x  = x0
    dt = 0.01_real64
 
    call minimize_fire(2, 2, x, dt, tol_grad, max_iter, &
      sphere_gradient_2d, write_iter_info, energy, ierr, mass, integrator)
 
    write(message(1), "(a,i6,a)") "FIRE (Sphere) converged in ", -ierr, " iterations"
    write(message(2), "(a,2(f8.4,a))")  " Minimum found at (", x(1), ", ", x(2), ")"
    write(message(3), "(a)")            " Analytical minimum is (0,0)."
    call messages_info(3)
    call messages_new_line()
 
 
    !=== Rastrigin ==========================================================
    call set_log_file("rastrigin.log")
    call io_rm(trim(current_log), namespace)
 
    x0 = [3.0_real64, 4.0_real64]
    x  = x0
    dt = 0.001_real64
    tol_grad = 1e-10_real64
 
    call minimize_fire(2, 2, x, dt, tol_grad, max_iter, &
      rastrigin_gradient_2d, write_iter_info, energy, ierr, mass, integrator)
 
    write(message(1), "(a,i6,a)") "FIRE (Rastrigin) converged in ", -ierr, " iterations"
    write(message(2), "(a,2(f8.4,a))")  " Minimum found at (", x(1), ", ", x(2), ")"
    write(message(3), "(a)")            " Analytical minimum is (0,0)."
    call messages_info(3)
 
    deallocate(current_log)
 
    pop_sub(test_optimizers)
  end subroutine test_optimizers
 
  !---------------------------------------------------------------------------
  subroutine rosenbrock_gradient_2d(n, x, val, getgrad, grad) !< get the new gradients given a new x
    integer,      intent(in)    :: n
    real(real64), intent(in)    :: x(n)
    real(real64), intent(inout) :: val
    integer,      intent(in)    :: getgrad
    real(real64), intent(inout) :: grad(n)
 
    real(real64), parameter :: a = m_one
    real(real64), parameter :: b = 100.0_real64
 
    push_sub(rosenbrock_gradient_2d)
 
    ! Only the 2D version is implemented here
    assert(n==2)
 
    val = (a-x(1))**2 +  b * (x(2) -x(1)**2)**2
 
    grad(1) = -m_two * (a-x(1)) - m_four * b * (x(2) -x(1)**2) * x(1)
    grad(2) =  m_two * b * (x(2) -x(1)**2)
 
    pop_sub(rosenbrock_gradient_2d)
  end subroutine rosenbrock_gradient_2d
 
 
  !---------------------------------------------------------------------------
  subroutine sphere_gradient_2d(n, x, val, getgrad, grad)
    integer, intent(in)          :: n
    real(real64), intent(in)     :: x(n)
    real(real64), intent(inout)  :: val
    integer, intent(in)          :: getgrad
    real(real64), intent(inout)  :: grad(n)
 
    push_sub(sphere_gradient_2d)
    assert(n==2)
 
    val     = dot_product(x, x)
    grad = m_two * x
 
    pop_sub(sphere_gradient_2d)
  end subroutine sphere_gradient_2d
 
 
  !---------------------------------------------------------------------------
  subroutine rastrigin_gradient_2d(n, x, val, getgrad, grad)
    integer, intent(in)          :: n
    real(real64), intent(in)     :: x(n)
    real(real64), intent(inout)  :: val
    integer, intent(in)          :: getgrad
    real(real64), intent(inout)  :: grad(n)
 
    real(real64), parameter :: a  = 10.0_real64
 
    push_sub(rastrigin_gradient_2d)
    assert(n==2)
 
    val = m_two*a + &
      (x(1)**2 - a*cos(m_two * m_pi * x(1))) + &
      (x(2)**2 - a*cos(m_two * m_pi * x(2)))
 
    grad = m_two*(x + m_pi * a * sin(m_two * m_pi * x))
 
    pop_sub(rastrigin_gradient_2d)
  end subroutine rastrigin_gradient_2d
 
  !---------------------------------------------------------------------------
  subroutine write_iter_info(iter, n, val, maxdr, maxgrad, x) !< output for each iteration step
    integer, intent(in) :: iter
    integer, intent(in) :: n
    real(real64), intent(in) :: val
    real(real64), intent(in) :: maxdr
    real(real64), intent(in) :: maxgrad
    real(real64), intent(in) :: x(n)
 
    integer :: iunit
 
    if (mpi_world%is_root()) then
      iunit = io_open(current_log, global_namespace, action = 'write', position = 'append')
      if (iter == 1) then
        write(iunit, '(a10, 3a20)') '#     iter','x','y', 'max_force'
      end if
      write(iunit, '(i10,3f20.10)') iter, x(1), x(2), maxgrad
      call io_close(iunit)
    end if
 
  end subroutine write_iter_info
 
  !---------------------------------------------------------------------------
  subroutine set_log_file(name)
    character(len=*), intent(in) :: name
    if (allocated(current_log)) deallocate(current_log)
    allocate(character(len=len_trim(name)) :: current_log)
    current_log = trim(name)
  end subroutine set_log_file
 
end module minimizer_tests_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: