Write algorithm for general usage

Help
1

Hi everyone, I am wondering whether it is possible to create an “interface” or “pointer” on the objective function that the algorithm wants to minimize.

The following is my implementation of golden section search. I am writing a dynamic programming problem, which requires the following algorithm to be given all expectation, idx, xdval, xuval four variables inside a for loop, and then using this implementation to compute the minimum (or maximum).

The part that I want to delegate for a user to modify is the interface block and somehow inside the gss_1d function I can somehow point or refer to whatever user-defined and calculate the value of the objective function.

module goldenSectionSearch

    ! ----------------------------------------------------- !
    ! modify the interface based on your objective function !
    ! (adjust numbers of arguments, intent and optional)    !
    ! and also whenever the corresponding gss subroutine    !
    ! calls the Obj procedure                               !
    ! ----------------------------------------------------- !

    use iso_Fortran_env, only: rk => real64, ik => int32
    implicit none
    private
    public :: gss

    real(rk), parameter :: goldenRatio = (3.0_rk - dsqrt(5.0_rk)) / 2.0_rk

    interface gss
        module procedure gss_1d
    end interface

    interface
        function gssObj(x, expectation, idx, xdval, xuval) result(f)
            import :: rk, ik
            ! class(*), intent(in), optional :: func_data
            integer(ik), intent(in) :: idx
            real(rk), dimension(:, :, :), intent(in) :: expectation
            real(rk), intent(in), optional :: xdval, xuval
            real(rk), intent(in) :: x
            real(rk) :: f
        end function gssObj
    end interface


contains

    subroutine gss_1d(fout, xout, Obj, LB, UB, expectation, idx, xdval, xuval, tol, maxiter, findmax, show)
        procedure(gssObj) :: Obj
        real(rk), intent(out) :: fout, xout
        real(rk), intent(in) :: LB, UB
        integer(ik), intent(in) :: idx
        real(rk), dimension(:, :, :), intent(in):: expectation
        real(rk), intent(in), optional :: xdval, xuval
        real(rk), intent(in), optional :: tol
        integer(ik), intent(in), optional :: maxiter
        logical, intent(in), optional :: findmax
        logical, intent(in), optional :: show
        ! class(*), intent(in), optional :: func_data
        character(len=128) :: tmp
        integer(ik) :: iter, gssMaxIter
        logical :: isFindMax, isShow
        real(rk) :: a, b, c, d, z, fval, fc, fd, gssTol, gssDist


        ! ------- !
        ! default !
        ! ------- !
        iter = 0_ik
        gssTol = 1D-12
        gssMaxIter = 500_ik
        isFindMax = .true.
        isShow = .false.
        if (present(tol)) gssTol = tol
        if (present(maxiter)) gssMaxIter = maxiter
        if (present(findmax)) isFindMax = findmax
        if (present(show)) isShow = show

        ! --------------------- !
        ! golden section search !
        ! --------------------- !

        gssDist = 2.0_rk*gssTol

        a = LB
        b = UB
        c = a + goldenRatio*(b-a);
        d = a + (1.0_rk-goldenRatio)*(b-a);

        fval = Obj(c, expectation, idx, xdval, xuval); fc = merge(-fval, fval, isFindMax);
        fval = Obj(d, expectation, idx, xdval, xuval); fd = merge(-fval, fval, isFindMax);

        do while (gssDist > gssTol .and. iter <= gssMaxIter)

            iter = iter + 1

            if (fc >= fd) then
                z = c + (1.0_rk-goldenRatio)*(b-c)
                ! case 1 [a c d b] <--- [c d z b]
                a = c
                c = d
                fc = fd
                d = z
                fval = Obj(d, expectation, idx, xdval, xuval); fd = merge(-fval, fval, isFindMax);
            else
                z = a + goldenRatio*(d-a)
                ! case 2 [a c d b] <--- [a z c d]
                b = d
                d = c
                fd = fc
                c = z
                fval = Obj(c, expectation, idx, xdval, xuval); fc = merge(-fval, fval, isFindMax);
            endif

            gssDist = b - a

            if (iter == 1 .and. isShow) then
                write(*, *) ''
                write(*, '(a4, 5(a20))') 'iter', '||b-a||', 'a', 'c', 'd', 'b'
                write(tmp, '(a4, 5(a20))') 'iter', '||b-a||', 'a', 'c', 'd', 'b'
                write(*, '(a)') repeat('=', len_trim(tmp))
            endif
            if (isShow) write (*, '(i4, 5(ES20.6))') iter, gssDist, a, c, d, b

            if (iter > gssMaxIter) write(*, *) "Warning: Max Iter reached."

        enddo

        if (z < LB) then
            z = LB
            fout = Obj(z, expectation, idx, xdval, xuval);
            xout = z
        elseif (z > UB) then
            z = UB
            fout = Obj(z, expectation, idx, xdval, xuval);
            xout = z
        else
            fout = fval
            xout = z
        endif


    end subroutine gss_1d

end module goldenSectionSearch
2

I don’t get your point… You would like it be possible for the user to provide an objective function with any interface, without having to modify your module? Am I correct?

3

The problem with just passing a procedure is that you have to hard code every parameter inside the objective function because Fortran doesn’t let you create function closures (yet?). That is often too restrictive. I would use an abstract type like this:

module objective_function_mod
    use iso_Fortran_env, only: dp => real64
    implicit none

    private
    public objective_function_t

    type, abstract :: objective_function_t
    contains
        procedure(eval_i), deferred eval
    end type

    abstract interface
        function eval_i(this, x) result(f)
            import :: dp, objective_function_t
            class(objective_function_t), intent(in) :: this
            real(dp), intent(in) :: x
            real(dp) :: f
        end function
    end interface
end module

Then you can implement a concrete objective function like this:

module quadratic_function_mod
    use iso_Fortran_env, only: dp => real64
    use objective_function_mod, only: objective_function_t
    implicit none

    private
    public quadratic_function_t

    type, extends(objective_function_t) :: quadratic_function_t
        real(rk) :: a, b, c
    contains
        procedure :: eval
    end type

contains

    function eval(this, x) result(f)
        class(objective_function_t), intent(in) :: this
        real(dp), intent(in) :: x
        real(dp) :: f

        f = this%a * x**2 + this%b * x + this%c
    end function

end module

Obviously you need to adapt the eval function in my example to your specific needs.

1 Like
4

You could do this by defining a module with a contained routine. Sketch:

module my_objective
     implicit none

    ... any data you need
    integer :: r, s, t ! Additional parameters

contains
subroutine find_minimum( param1, param2, ... )
      call golden_seection( xmin,xmax, my_obj )
contains
real function my_obj( x )
     my_obj = func( param1, param2, x )|
end function my_obj
end subroutine find_minimum

! The actual function ...
real func( p, q, x )
     func = complicated_calculation( p, q, r, s, t, x )
end function func
end module my_objective

This way you have access to any “global” data (r,s t) and arguments that you need (param1, param2) without the interface to the minimization routine restricting you.

As others have shown, there are other solutions as well.

1 Like
5

@fish830911 You may also consider reverse communication techniques to write a kind of generic algorithm.

6

@PierU Yes that’s what I mean, sorry that I am not clear on what is my objective. I just don’t want to modify the objective inside the algorithm.

@plevold May I ask in your example, is the eval function corresponds to the objective function, or the algorithm that I am trying to implement?

@Arjen May I ask in this sketch, I need to define additional parameters r, s, t inside the module my_objective?

8

The objective function, so you’d need to adapt it to fit your gssObj interface. You can then use the abstract type like so:

subroutine gss_1d(objective_function) ! Add other arguments as needed
    class(objective_function_t), intent(in) :: objective_function

    real(dp) :: f
    real(dp) :: x
 
    ! ...
    ! Evaluate objective function
    f = objective_function%eval(x)
    ! ...
end subroutine
1 Like
9

No, that was merely an example to show the possibilities.

10

Just curious, will this OOP method works with openmp if the additional parameters (I guess is a, b, c in your example) is private inside the do loop?

11

I see no reason why it shouldn’t. Multithreading like OpenMP with shared state (which a, b and c will be in this case) can become problematic when that state is mutated. We guard against this by requiring the this argument to the eval function to be intent(in). If possible you can try to make the function pure as well. This will guard against mutating any global state during evaluation of the objective function.

1 Like
12

I fixed a few minor issues in your example and it works… I then tried to put the content of the quadratic_function_mod module directly in a main program like this:

program foo
    use iso_Fortran_env, only: dp => real64
    use objective_function_mod

    type, extends(objective_function_t) :: quadratic_function_t
        real(dp) :: a, b, c
    contains
        procedure :: eval
    end type

    continue ! does nothing at the moment

contains

    function eval(this, x) result(f)
        class(quadratic_function_t), intent(in) :: this
        real(dp), intent(in) :: x
        real(dp) :: f

        f = this%a * x**2 + this%b * x + this%c
    end function

end program

and I get the following errors:

   57 |         procedure :: eval
      |                 1
Error: 'eval' must be a module procedure or an external procedure with an explicit interface at (1)

   54 |     type, extends(objective_function_t) :: quadratic_function_t
      |                                                               1
Error: Derived-type 'quadratic_function_t' declared at (1) must be ABSTRACT because 'eval' is DEFERRED and not overridden

It looks like the actual code of the eval() has to be in a module. Why not, but what is the technical reason behind that ?

1 Like
13

I have came up with another solution without type bound procedure (not better or whatever, just different). Rewriting the solution of @plevold (an after a review from @ivanpribec to simplify it):

The general algorithm

module objective_function_mod
    implicit none

    private
    public optimize

    abstract interface
        real function userfunc_i(p, x)
            class(*), intent(in) :: p
            real, intent(in) :: x
        end function
    end interface

contains
   
   subroutine optimize(userfunc,p,x0,x)
      procedure(userfunc_i) :: userfunc
      class(*), intent(in) :: p
      real, intent(in) :: x0
      real, intent(out) :: x
      ! ...
      ! some code with calls to the objective function :
      ! f = userfunc(p,x)
      ! ...  
   end subroutine

end module

The specific function and its parameters:

module quadratic_function_mod
    use objective_function_mod
    implicit none

    private
    public myparams_t, myfunc

    type myparams_t
        real :: a, b, c
    end type

contains

    real function myfunc(p, x) result(f)
        class(*), intent(in) :: p
        real, intent(in) :: x

        select type(p)
            type is (myparams_t)     
                f = p%a * x**2 + p%b * x + p%c
       end select
   end function

end module

And how to use:

...
use objective_function_mod
use quadratic_function_mod
...
x0 = ...
call optimize(myfunc,                    &
              myparams_t(3.0,-2.0,-1.0), &   ! a, b, c
              x0,                        &
              x                          )
...
14

Following a remark from @ivanpribec, a generic class(*) can be used to further simplify my solution above (edited)

15

Fortran functions don’t have parameters. They have arguments.

17

Sounds good, but I was talking about the conceptual objective function.