Merge pull request #26 from jacobwilliams/25-bounds

added a simple bounds check mode

Author: jacobwilliams

.github/workflows/CI.yml

@@ -10,8 +10,8 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest]
-        gcc_v: [10] # Version of GFortran we want to use.
-        python-version: [3.9]
+        gcc_v: [14] # Version of GFortran we want to use.
+        python-version: [3.12]
     env:
       FC: gfortran-${{ matrix.gcc_v }}
       GCC_V: ${{ matrix.gcc_v }}
@@ -31,7 +31,7 @@ jobs:
       uses: ts-graphviz/setup-graphviz@v1
 
     - name: Setup Fortran Package Manager
-      uses: fortran-lang/setup-fpm@v5
+      uses: fortran-lang/setup-fpm@v7
       with:
         github-token: ${{ secrets.GITHUB_TOKEN }}
 
@@ -53,6 +53,7 @@ jobs:
            --install /usr/bin/gcc gcc /usr/bin/gcc-${{ matrix.gcc_v }} 100 \
            --slave /usr/bin/gfortran gfortran /usr/bin/gfortran-${{ matrix.gcc_v }} \
            --slave /usr/bin/gcov gcov /usr/bin/gcov-${{ matrix.gcc_v }}
+        sudo apt-get install libblas-dev liblapack-dev
 
     # - name: Compile
     #   run: fpm build --profile release
@@ -80,7 +81,8 @@ jobs:
 
     - name: Deploy Documentation
       if: github.ref == 'refs/heads/master'
-      uses: JamesIves/github-pages-deploy-action@v4.4.1
+      uses: JamesIves/github-pages-deploy-action@v4.7.3
       with:
         branch: gh-pages # The branch the action should deploy to.
         folder: doc  # The folder the action should deploy.
+        single-commit: true

README.md

@@ -31,6 +31,9 @@ A work in progress.
      * backtracking linesearch method
      * exact linesearch method using [fmin](https://github.com/jacobwilliams/fmin) minimizer
      * evaluate function at specified fixed points
+  * Has two options for variable bounds (`xlow<=x<=xupp`):
+     * Ignore bounds
+     * Crude method: manually adjust `x` vector at each function evaluation so that `x = min(max(x,xlow),xupp)`.
 
 ### Compiling
 

src/nlesolver_module.F90

@@ -108,6 +108,15 @@ module nlesolver_module
         character(len=:),allocatable :: message         !! latest status message
         integer     :: istat            = -999          !! latest status message
 
+        integer :: bounds_mode = 0  !! how to handle the `x` variable bounds:
+                                    !!
+                                    !!  * 0 = ignore bounds.
+                                    !!  * 1 = use bounds (if specified) by adjusting the `x` vector
+                                    !!    at each function evaluation so that each individual `x`
+                                    !!    component is within its bounds.
+        real(wp),dimension(:),allocatable :: xlow !! lower bounds for `x` (size is `n`). only used if `bounds_mode>0`.
+        real(wp),dimension(:),allocatable :: xupp !! upper bounds for `x` (size is `n`). only used if `bounds_mode>0`.
+
         procedure(func_func),pointer       :: func             => null() !! user-supplied routine to compute the function
         procedure(export_func),pointer     :: export_iteration => null() !! user-supplied routine to export iterations
         procedure(wait_func),pointer       :: user_input_check => null() !! user-supplied routine to enable user to stop iterations
@@ -167,6 +176,7 @@ module nlesolver_module
         procedure,public :: status     => get_status
 
         procedure :: set_status
+        procedure :: adjust_x_for_bounds
 
         end type nlesolver_type
     !*********************************************************
@@ -270,27 +280,14 @@ subroutine set_status(me,istat,string,i,r)
     integer,intent(in),optional         :: i       !! an integer value to append
     real(wp),intent(in),optional        :: r       !! a real value to append
 
-    character(len=256) :: numstr !! for number fo string conversion
     character(len=:),allocatable :: message !! the full message to log
     integer :: iostat !! write `iostat` code
 
     message = trim(string)
-    if (present(i)) then
-        numstr = ''
-        write(numstr,fmt=*,iostat=iostat) i
-        if (iostat/=0) numstr = '****'
-        message = message//' '//trim(adjustl(numstr))
-    end if
-    if (present(r)) then
-        numstr = ''
-        write(numstr,fmt=*,iostat=iostat) r
-        if (iostat/=0) numstr = '****'
-        message = message//' '//trim(adjustl(numstr))
-    end if
+    if (present(i)) message = message//' '//int2str(i)
+    if (present(r)) message = message//' '//real2str(r)
 
-    if (me%verbose) then
-        write(me%iunit,'(A)',iostat=iostat) message
-    end if
+    if (me%verbose)  write(me%iunit,'(A)',iostat=iostat) message
 
     ! store in the class:
     me%istat = istat
@@ -299,6 +296,42 @@ subroutine set_status(me,istat,string,i,r)
     end subroutine set_status
 !*****************************************************************************************
 
+!*****************************************************************************************
+!>
+!  Convert an integer to a string. Works for up to 256 digits.
+
+    function int2str(i) result(s)
+        integer, intent(in) :: i !! integer to convert
+        character(len=:),allocatable :: s !! string result
+        character(len=256) :: tmp !! temp string
+        integer :: iostat !! write `iostat` code
+        write(tmp,fmt=*,iostat=iostat) i
+        if (iostat/=0) then
+            s = '****'
+        else
+            s = trim(adjustl(tmp))
+        end if
+    end function int2str
+!*****************************************************************************************
+
+!*****************************************************************************************
+!>
+!  Convert a real to a string. Works for up to 256 digits.
+
+    function real2str(r) result(s)
+        real(wp), intent(in) :: r !! real to convert
+        character(len=:),allocatable :: s !! string result
+        character(len=256) :: tmp !! temp string
+        integer :: iostat !! write `iostat` code
+        write(tmp,fmt=*,iostat=iostat) r
+        if (iostat/=0) then
+            s = '****'
+        else
+            s = trim(adjustl(tmp))
+        end if
+    end function real2str
+!*****************************************************************************************
+
 !*****************************************************************************************
 !>
 !  Return the status code and message from the [[nlesolver_type]] class.
@@ -320,6 +353,7 @@ end subroutine set_status
 !  * -13  -- Error: backtracking linesearch tau must be in range (0, 1)
 !  * -14  -- Error: must specify grad_sparse, irow, and icol for sparsity_mode > 1
 !  * -15  -- Error: irow and icol must be the same length
+!  * -16  -- Error: xlow > xupp
 !  * -999 -- Error: class has not been initialized
 !  * 0    -- Class successfully initialized in [[nlesolver_type:initialize]]
 !  * 1    -- Required accuracy achieved
@@ -361,7 +395,8 @@ subroutine initialize_nlesolver_variables(me,&
                     verbose,iunit,n_uphill_max,n_intervals,&
                     sparsity_mode,irow,icol,&
                     atol,btol,conlim,damp,itnlim,nout,&
-                    lusol_method,localSize,custom_solver_sparse )
+                    lusol_method,localSize,custom_solver_sparse,&
+                    bounds_mode,xlow,xupp)
 
     implicit none
 
@@ -437,6 +472,16 @@ subroutine initialize_nlesolver_variables(me,&
                                              !! At most `min(m,n)` vectors will be allocated.
     procedure(sparse_solver_func),optional :: custom_solver_sparse !! for `sparsity_mode=5`, this is the
                                                                    !! user-provided linear solver.
+    integer,intent(in),optional :: bounds_mode  !! how to handle the `x` variable bounds:
+                                                !!
+                                                !!  * 0 = ignore bounds
+                                                !!  * 1 = use bounds (if specified) by adjusting the `x` vector
+                                                !!    at each step so that each individual `x` component is within
+                                                !!    the bounds
+    real(wp),dimension(n),intent(in),optional :: xlow  !! lower bounds for `x` (size is `n`). only used if `bounds_mode>0` and
+                                                       !! both `xlow` and `xupp` are specified.
+    real(wp),dimension(n),intent(in),optional :: xupp  !! upper bounds for `x` (size is `n`). only used if `bounds_mode>0` and
+                                                       !! both `xlow` and `xupp` are specified.
 
     logical :: status_ok !! true if there were no errors
 
@@ -453,6 +498,19 @@ subroutine initialize_nlesolver_variables(me,&
 
     !optional:
 
+    if (present(bounds_mode) .and. present(xlow) .and. present(xupp)) then
+        if (any(xlow>xupp)) then  ! check for consistency
+            status_ok = .false.
+            call me%set_status(istat = -16, string = 'Error: xlow > xupp')
+            return
+        end if
+        me%bounds_mode = bounds_mode
+        me%xupp = xupp
+        me%xlow = xlow
+    else
+        me%bounds_mode = 0  ! default
+    end if
+
     if (present(step_mode)) then
         select case (step_mode)
         case(1) ! = use the specified `alpha` (0,1]
@@ -672,6 +730,7 @@ subroutine nlesolver_solver(me,x)
     end if
 
     ! evaluate the function:
+    call me%adjust_x_for_bounds(x) ! if the guess is out of bounds it may also be adjusted first.
     call me%func(x,fvec)
     f = norm2(fvec)
 
@@ -923,6 +982,35 @@ subroutine nlesolver_solver(me,x)
     end subroutine nlesolver_solver
 !*****************************************************************************************
 
+!*****************************************************************************************
+!>
+!  if necessary, adjust the `x` vector to be within the bounds.
+
+    subroutine adjust_x_for_bounds(me,x)
+
+    implicit none
+
+    class(nlesolver_type),intent(inout) :: me
+    real(wp),dimension(me%n),intent(inout) :: x  !! the `x` vector to adjust
+
+    integer :: i !! counter
+
+    if (me%bounds_mode==1) then
+        ! x = min(max(x,me%xlow),me%xupp)
+        do i = 1, me%n
+            if (x(i)<me%xlow(i)) then
+                x(i) = me%xlow(i)
+                if (me%verbose) write(me%iunit, '(A)') 'x('//int2str(i)//') < xlow(i) : adjusting to lower bound'
+            else if (x(i)>me%xupp(i)) then
+                x(i) = me%xupp(i)
+                if (me%verbose) write(me%iunit, '(A)') 'x('//int2str(i)//') > xupp(i) : adjusting to upper bound'
+            end if
+        end do
+    end if
+
+    end subroutine adjust_x_for_bounds
+!*****************************************************************************************
+
 !*****************************************************************************************
 !>
 !   Destructor
@@ -1042,6 +1130,7 @@ subroutine simple_step(me,xold,p,x,f,fvec,fjac,fjac_sparse)
     real(wp),dimension(:),intent(in),optional :: fjac_sparse !! jacobian matrix [sparse]
 
     x = xold + p * me%alpha
+    call me%adjust_x_for_bounds(x)
 
     !evaluate the function at the new point:
     call me%func(x,fvec)
@@ -1115,6 +1204,7 @@ subroutine backtracking_linesearch(me,xold,p,x,f,fvec,fjac,fjac_sparse)
     do
 
         xtmp = xold + p * alpha
+        call me%adjust_x_for_bounds(xtmp)
         call me%func(xtmp,fvectmp)
         ftmp = norm2(fvectmp)
 
@@ -1183,6 +1273,7 @@ subroutine exact_linesearch(me,xold,p,x,f,fvec,fjac,fjac_sparse)
     if (me%verbose) write(me%iunit,'(1P,*(A,1X,E16.6))') '        alpha_min = ', alpha_min
 
     x = xold + p * alpha_min
+    call me%adjust_x_for_bounds(x)
     if (all(x==xnew)) then
         ! already computed in the func
     else
@@ -1198,6 +1289,7 @@ real(wp) function func_for_fmin(alpha)
     real(wp),intent(in) :: alpha !! indep variable
 
     xnew = xold + p * alpha
+    call me%adjust_x_for_bounds(xnew)
     call me%func(xnew,fvec)
     func_for_fmin = norm2(fvec) ! return result
 
@@ -1260,6 +1352,7 @@ subroutine fixed_point_linesearch(me,xold,p,x,f,fvec,fjac,fjac_sparse)
     do i = 1, n_points
 
         x_tmp = xold + p * alphas_to_try(i)
+        call me%adjust_x_for_bounds(x_tmp)
 
         ! evaluate the function at tthis point:
         call me%func(x_tmp,fvec_tmp)

test/nlesolver_test_1.f90

@@ -98,7 +98,10 @@ program nlesolver_test_1
                                 export_iteration = export,&
                                 n_intervals = n_intervals, &
                                 fmin_tol = fmin_tol, &
-                                verbose = verbose)
+                                verbose = verbose, &
+                                bounds_mode = 1, &
+                                xlow = [0.0_wp, -5.0_wp], &
+                                xupp = [1.0_wp, 0.0_wp])
         call solver%status(istat, message)
         write(*,'(I3,1X,A)') istat, message
         if (istat /= 0) error stop