(Fortran) Using Procedures as variables for Derived data type

[Heo__Jun-Yeong]

Hi, I'm working on optimizing my codes written in Fortran. 

In my working field, there is famous python package, 'SimPEG'.

(https://docs.simpeg.xyz/index.html  -- main homepage)

(https://github.com/simpeg/simpeg -- github)

The following code is possible in python with SimPEG

...

dmis_grav = data_misfit.L2DataMisfit(data=data_object_grav, simulation=simulation_grav)
dmis_mag = data_misfit.L2DataMisfit(data=data_object_mag, simulation=simulation_mag)

...

dmis = dmis_grav + dmis_mag

...

inv_prob = inverse_problem.BaseInvProblem(dmis, reg, opt)

...

inv = inversion.BaseInversion(inv_prob, directives_list)

In the upper code, Inversion is optimization process (i.e. Gauss-Newton), and iteratively calls 'dmis' to calculate data misfit for different models.

The main point is we can make new function (or class?) by simple statement:

dmis = dmis_grav + dmis_mag.

(Everytime we call 'dmis', the program calls 'dmis_grav' and 'dmis_mag', and add them)

So, I thinking about:

Is there any coding techniques which can make this possible in Fortran?

The simple code of my idea is:

module mdata
        implicit none
        type :: tdat
          real, dimension(2) :: d
        end type tdata
end module mdata

module mmodel
        implicit none
        type :: tmodel
          real, dimension(2) :: m
        end type tmodel
end module mmodel

module msimulation1
        use mdata
        use mmodel
        implicit none
        type :: tsimulation1
          type(tmodel) :: model
         contains
          procedure :: init_simulation1
          procedure :: simulation1
        end type tsimulation1
  contains
    subroutine init_simulation1(this, inmodel)
        class(tsimulation1) :: this
        type(tmodel), intent(in) :: inmodel

        this%model = inmodel
    end subroutine init_simulation1

    subroutine simulation1(this, outdata)
        class(tsimulation1), intent(in) :: this
        type(tdata), intent(out) :: outdata

        outdata%d(:) = this%model%m(:)*2.
    end subroutine simulation1
end module msimulation1

module msimulation2
        use mdata
        use mmodel
        implicit none
        type :: tsimulation2
          type(tmodel) :: model
         contains
          procedure :: init_simulation2
          procedure :: simulation2
        end type tsimulation2
  contains
    subroutine init_simulation2(this, inmodel)
        class(tsimulation2) :: this
        type(tmodel), intent(in) :: inmodel

        this%model = inmodel
    end subroutine init_simulation2

    subroutine simulation2(this, outdata)
        class(tsimulation1), intent(in) :: this
        type(tdata), intent(out) :: outdata

        outdata%d(:) = this%model%m(:)/2.
    end subroutine simulation2
end module msimulation2

for the initial setting, 

and need to make misfit module such as:

module mMisfit
        use mdata
        implicit none
        type :: tL2
          type(tdata) :: d
          real, external :: s
         contains
          procedure :: init_L2
          procedure :: L2DataMisfit
        end type tL2
  contains
    subroutine init_L2(this, indata, insimulation)
        class(tL2) :: this
        type(tdata), intent(in) :: indata
        real, external :: insimulation
        
        this%d = indata
        this%s = insimulation
    end subroutine init_L2

    real function L2DataMisfit(this)
        class(tL2) :: this
        type(tdata) :: tmpd
        real :: sdot

        call this%s(tmpd)

        L2DataMisfit = sdot(2,this%d%d-tmpd%d,1)
    end function L2DataMisfit
end module mMisfit

The upper module is impossible because of the

'real, external :: s' in 'type :: tL2'.

However, this kind of procedure-variable (?) or procedure-passing (?) is necessary to make the 'L2DataMisfit' function run for both 'simulation1' and 'simulation2' without changing module code, and only dependent on main program. 

Is there any good idea to make this possible in Fortran...?

[jimdempseyatthecove · Accepted solution]

Use a procedure pointer.

The example given in the link above is not complete.

The first example should have declared the interface to (or of) QUARK.

Note, the interface can be declared as an ABSTRACE INTERFACE should you not have declared a specific actual interface.

Jim Dempsey

[Heo__Jun-Yeong]

Thank you for your reply. 

I haven't known about 'ABSTRACT' or 'procedure pointer' in Fortran and I studied about them.

Below is the fixed version of my code

module mdata
        implicit none
        type :: tdat
          real, dimension(2) :: d
        end type tdata
end module mdata

module mmodel
        implicit none
        type :: tmodel
          real, dimension(2) :: m
        end type tmodel
end module mmodel

module msimulation
        use mdata
        use mmodel
        implicit none
        type :: tsimulation
          type(tmodel) :: model
         contains
          procedure(init_simulation), deferred :: init
          procedure(simulation), deferred :: run
        end type tsimulation

        abstract interface
          subroutine init_simulation(this, inmodel)
            import tsimulation, tmodel
            class(tsimulation) :: this
            type(tmodel) :: inmodel
          end subroutine init_simulation

          subroutine simulation(this, outdata)
            import tsimulation, tdata
            class(tsimulation) :: this
            type(tdata) :: outdata
          end subroutine simulation
        end interface
end module msimulation

module msimulation1
        use msimulation
        implicit none
        type, extends(simulation) :: tsimulation1
         contains
          procedure :: init => init_simulation1
          procedure :: run => simulation1
        end type tsimulation1
  contains
    subroutine init_simulation1(this, inmodel)
        class(tsimulation1) :: this
        type(tmodel), intent(in) :: inmodel

        this%model = inmodel
    end subroutine init_simulation1

    subroutine simulation1(this, outdata)
        class(tsimulation1), intent(in) :: this
        type(tdata), intent(out) :: outdata

        outdata%d(:) = this%model%m(:)*2.
    end subroutine simulation1
end module msimulation1

module msimulation2
        use msimulation
        implicit none
        type, extends(simulation) :: tsimulation2
         contains
          procedure :: init => init_simulation2
          procedure :: run => simulation2
        end type tsimulation2
  contains
    subroutine init_simulation2(this, inmodel)
        class(tsimulation2) :: this
        type(tmodel), intent(in) :: inmodel

        this%model = inmodel
    end subroutine init_simulation2

    subroutine simulation2(this, outdata)
        class(tsimulation2), intent(in) :: this
        type(tdata), intent(out) :: outdata

        outdata%d(:) = this%model%m(:)/2.
    end subroutine simulation2
end module msimulation2

and misfit is

module mMisfit
        use mdata
        use msimulation
        implicit none
        type :: tL2
          type(tdata) :: d
          class(tsimulation), pointer :: s => null()
         contains
          procedure :: init_L2
          procedure :: L2DataMisfit
        end type tL2
  contains
    subroutine init_L2(this, indata, insimulation)
        class(tL2) :: this
        type(tdata), intent(in) :: indata
        class(tsimulation), target :: insimulation
        
        this%d = indata
        this%s => insimulation
    end subroutine init_L2

    real function L2DataMisfit(this)
        class(tL2) :: this
        type(tdata) :: tmpd
        real :: sdot

        call this%s%run(tmpd)

        L2DataMisfit = sdot(2,this%d%d-tmpd%d,1)
    end function L2DataMisfit
end module mMisfit

I used pointer for derived data type 'tL2' because in my future work, the 'simulation' class does not need to be copied (and it might require large memories).

And it works fine as:

program Main
        use mdata
        use mmodel
        use mMisfit
        use msimulation1
        implicit none
        type(tdata) :: d0
        type(tmodel) :: m0
        type(tsimulation1) :: s1
        type(tL2) :: dmis

        d0%d(:) = [10., 10.]
        m0%m(:) = [2., 3.]
        call s1%init(m0)
        call dmis%init_L2(d0, s1)
        print*, dmis%L2DataMisfit()
end program Main

Now, to the original question, 

In python, we can make new class(?) simply by:

dmis_grav = data_misfit.L2DataMisfit(data=data_object_grav, simulation=simulation_grav)
dmis_mag = data_misfit.L2DataMisfit(data=data_object_mag, simulation=simulation_mag)

dmis = dmis_grav + dmis_mag

which can be used like:

inv_prob1 = inverse_problem.BaseInvProblem(dmis_mag, reg, opt)
inv_prob2 = inverse_problem.BaseInvProblem(dmis_grav, reg, opt)
inv_prob3 = inverse_problem.BaseInvProblem(dmis, reg, opt)

I think it is impossible in fortran...

The only way is defining new derived data type by extension something like:

type, abstract :: tmisfit
 contains
  procedure(init_misfit), deferred :: init
  procedure(DataMisfit), deferred :: misfit
end type tmisfit

type, extends(tmisfit) :: tL2
  type(tdata) :: d
  class(tsimulation), pointer :: s => null()
 contains
  procedure :: init => init_L2
  procedure :: misfit => L2DataMisfit
end type tL2

type, extends(tmisfit) :: tL2plusL2
  type(tL2) :: L2_1, L2_2
 contains
  procedure :: init => init_L2plusL2
  procedure :: misfit => L2plusL2DataMisfit
end type tL2plusL2

Do you have any better ideas...?

[jimdempseyatthecove]

Perhaps you can look at using a generic interface.

module YourModule
    type tObject_grav
        real :: zzz
    end type tObject_grav
    
    type tObject_mag
        real :: zzz
    end type tObject_mag
    
    type tObject_dmis
        real :: zzz
    end type tObject_dmis
    
    ! Provide a generic name to the interface
    interface BaseInvProblem
        ! Specify the calling signatured procedures
        procedure :: BaseInvProblem_grav, BaseInvProblem_mag, BaseInvProblem_dmis
    end interface BaseInvProblem
contains
    ! Then called procedures with unique calling signatures
    function BaseInvProblem_grav(dmis, reg, opt) result(ret)
        type(tObject_grav) :: dmis
        real :: reg
        integer :: opt
        ! (code)
        ret = 0 ! return value
    end function BaseInvProblem_grav
    
    function BaseInvProblem_mag(dmis, reg, opt) result(ret)
        type(tObject_mag) :: dmis
        real :: reg
        integer :: opt
        ! (code)
        ret = 0 ! return value
    end function BaseInvProblem_mag
    
    function BaseInvProblem_dmis(dmis, reg, opt) result(ret)
        type(tObject_dmis) :: dmis
        real :: reg
        integer :: opt
        ! (code)
        ret = 0 ! return value
    end function BaseInvProblem_dmis
    
end module YourModule
    
program Console20
    use YourModule
    implicit none
    type(tObject_grav), allocatable :: dmis_mag(:)
    type(tObject_mag),  allocatable :: dmis_grav(:)
    type(tObject_dmis), allocatable :: dmis(:)
    real :: reg
    integer :: I, what, opt
!    ...
    do I=1, size(dmis_grav)
        ! ...
        what = BaseInvProblem(dmis_grav(I), reg, opt) + BaseInvProblem(dmis_mag(I), reg, opt) + BaseInvProblem(dmis(i), reg, opt)
        ! ...
    end do
    ! ...
    end program Console20

Something along the above sketch.

There are other ways to do this too.

Jim Dempsey

[jimdempseyatthecove]

I caution you that for large numbers of like data with high iteration counts (e.g. simulation models), that using OOP-style programming suffers with performance. With OOP-style (Array of Structures) the code may be traversing levels of code to get to the necessary operations on scalar variables. Using "traditional" Fortran (Structure of Arrays), the code generally has (relatively) immediate access to vectors of data suitable for SIMD operations (Single Instructions on Multiple Data). SoA typically has much better cache utilization as well.

Jim Dempsey