Using FORTRAN Intel MKL to build curve_fit() available in scipy.optimize

[MohammadShafieipour]

I would like to implement a curve fitting procedure using Intel MKL in FORTRAN. To make it clear, I coded up an equivalent python code as follows:

###start of python code for curve fitting###

from scipy.optimize import curve_fit

#The function to be curve fitted
def func(x,c1,c2,c3):
   return c1 + c2 * x ** (c3)

#x datapoints
xdat=[1,2,3,4,5,6,7,8,9,10]

#y datapoints
ydat=[4,13,28,49,76,109,148,193,244,301]

#initial guess
ig=[0.,5.0,1.3]

popt, pcov = curve_fit(func,xdat, ydat, ig)
print('popt=',popt)
print('pcov=',pcov)

###end of python code for curve fitting###

As seen above, I have 10 x and y data points (xdat,ydat) that I would like to fit according to func() and find the corresponding coefficient values (c1,c2,c3). We know that the solution to this problem is exactly c1=1, c2=3, and c3=2. Accordingly, the python code prints the output below:

popt= [1. 3. 2.]
pcov= [[ 2.65907915e-28 -3.31803210e-29 4.52198404e-30]
[-3.31803210e-29 6.08357690e-30 -8.63287134e-31]
[ 4.52198404e-30 -8.63287134e-31 1.23305020e-31]]

In order to implement this using MKL, I am looking at 2 examples in Intel directory: ".\Intel\MKL\10.0.2.019\examples\solver\source", namely:

• ex_nlsqp_f.f
• ex_nlsqp_bc_f.f

However, they both have the same number of "FUNCTION VARIABLES (N)" and "DIMENSION OF FUNCTION VALUE (M)". So I am not sure if these solutions are applicable to the problem I am trying to solve. Also, the available EXTENDET_POWELL subroutine, does not take x and coefficient values (e.g. c1,c2,c3) which makes it difficult to relate to the solution code up in python.

Could you please provide guidelines as to how I can modify ex_nlsqp_f.f or ex_nlsqp_bc_f.f in order to implement the curve fitting example shown above in FORTRAN using MKL. Or maybe there are better MKL examples for this purpose?

Note: In order to keep things simple, I did not introduce bounds in the above python code, but my final goal is to include lower and upper bounds as used in ex_nlsqp_bc_f.f.

Thank you in advance.

[ArpitaP_Intel]

Hi,

Thanks for reaching out to us.

We are working on it internally. We will get back to you.

Regards,

Arpita

[mecej4]

Here is the code to obtain the best coefficients, which is an adaptation of the  DTRNLSP example that comes with MKL. You will need to add code to compute and print the covariance matrix, if you want that. Similarly, if you want to specify bounds for the fit parameters, you will need to modify the code by calling dtrnlsp_bc, etc.

      module pdata
         implicit none
         integer, parameter :: m=10, n=3 !no. of data pairs, parameters to fit
         double precision, dimension(m) :: 
     +       xdat = [1,2,3,4,5,6,7,8,9,10],
     +       ydat = [4,13,28,49,76,109,148,193,244,301]
      end module pdata
      program nlfit
        use pdata, only : m, n
        use mkl_rci
        implicit none
        external            func
        integer             rci_request, iter,itrn,st_cr,successful
        integer             iter1, iter2, i, j, info(6)
        double precision    x(n), eps(6), jac_eps, fvec(m), fjac(m, n)
        double precision    rs, r1, r2
        type(handle_tr)     handle
!             
        do i = 1, 6
            eps (i) = 1.d-5
        end do
        iter1 = 1000
        iter2 = 100
        rs = 100.d0
        jac_eps = 1.d-8
        x = [0d0,5d0,1.3d0]
        do i = 1, m
            fvec (i) = 0.d0
            do j = 1, n
                fjac (i, j) = 0.d0
            end do
        end do
        if (dtrnlsp_init (handle, n, m, x, eps, iter1, iter2, rs)
     &      .ne. tr_success) then
            print *, ' Error in dtrnlsp_init'
            call mkl_free_buffers
            stop 1
        end if
        if (dtrnlsp_check (handle, n, m, fjac, fvec, eps, 
     &     info) .ne. tr_success) then
            print *, ' Error in dtrnlspbc_init'
            call mkl_free_buffers
            stop 1
        else
            if( info(1) .ne. 0 .or. 
     &              info(2) .ne. 0 .or. 
     &              info(3) .ne. 0 .or. 
     &              info(4) .ne. 0 ) then
                print *, ' Input parameters are not valid'
                call mkl_free_buffers
                stop 1
            end if
        end if
        rci_request = 0
        successful = 0
        itrn = 0
        print *,' itrn  ||f||'
        do while (successful == 0)
            if (dtrnlsp_solve (handle, fvec, fjac, rci_request)
     &              .ne. tr_success) then
                print *, ' Error in dtrnlsp_solve'
                call mkl_free_buffers
                stop 1
            end if
            select case (rci_request)
            case (-1, -2, -3, -4, -5, -6)
                successful = 1
            case (1)
                itrn = itrn+1
                call func (m, n, x, fvec)
                print '(i3,es10.2)',itrn,norm2(fvec)
            case (2)
                if (djacobi (func, n, m, fjac, x, jac_eps)
     &                  .ne. tr_success) then
                    print *, ' Error in djacobi'
                    call mkl_free_buffers
                    stop 1
                end if
            endselect
        end do
        if (dtrnlsp_get (handle, iter, st_cr, r1, r2)
     &          .ne. tr_success) then
            print *, ' Error in dtrnlsp_get'
            call mkl_free_buffers
            stop 1
        end if
        if (dtrnlsp_delete (handle) .ne. tr_success) then
            print *, ' Error in dtrnlsp_delete'
            call mkl_free_buffers
            stop 1
        end if

        call mkl_free_buffers
        if (r2 .lt. 1.d-5) then
            PRINT *, 'Best Fit Coefficients'
            print '(3ES13.5)',x(1:n)
            STOP 0
        else
            PRINT *, ' Dtrnlsp failed'
            stop 1
        end if
      end program nlfit

      subroutine func (m, n, x, f)
        use pdata, only : xdat,ydat
        implicit none
        integer m, n
        double precision x (*), f (*)
        integer i

        do i = 1, m
           f(i) = x(1)+x(2)*xdat(i)**x(3)-ydat(i)
        end do

      end subroutine func

[MohammadShafieipour]

Thank you very much for providing the code. I tried to compile it but there was a complain about (use mkl_rci). In the MKL example files, I have seen the use of (INCLUDE "mkl_rci.fi") which I can compile on my machine. So I am sure there is a way to prepare mkl_rci as a module and compile your code as is. Could you please let me know what is the best way to do that. 

[mecej4]

In the MKL include directory, there is a file called mkl_rci.f90. Compiling that file will generate the module files mkl_rci.mod and mkl_rci_type.mod.

You can copy those files to the MKL include\<arch> directory, where <arch> is either "ia32" or "intel64", as appropriate.

The net effect is the same as adding the line 

   INCLUDE "mkl_rci.fi"

to your source file, but avoids needless recompilation of the module.

[MohammadShafieipour]

I am looking at the following directory

C:\Program Files (x86)\Intel\MKL\10.0.2.019\include

I do see the following files in it

• mkl_rci.fi
• mkl_rci.h

but I do not find mkl_rci.f90 in it. Should I be looking for it elsewhere?

[mecej4]

My goodness, you are using a version of MKL that is at least a decade old! Try the following.

Create a short file, in fixed form Fortran (give it any convenient name), containing the following three lines, and compile it. 

c
      module mkl_rci
      include 'mkl_rci.fi'
      end module

 You may need to modify other parts of my code to make it work with that old version of MKL. For instance, the type handle_tr probably did not exist in that version, so you may need to replace it with INTEGER or INTEGER*8, depending on whether you are compiling for 32-bit or 64-bit. You may have to modify the ex_nlsqp_f.f version that you have, using my source code only as a guideline.

[MohammadShafieipour · Accepted solution]

Below are my findings. Based on mecej4 solution, I prepared the following f90 code. If you have already installed "Intel Parallel Studio" and included MKL in Visual Studio, it should compile as is and it produces the correct solution as follows:

itrn ||f||
1 3.02E+02
2 1.93E+04
3 4.27E+02
4 2.50E+02
5 3.74E+01
6 9.43E+00
7 7.82E+00
8 5.03E-01
9 6.13E-02
10 1.42E-06
Best Fit Coefficients
1.00000E+00 3.00000E+00 2.00000E+00

Note that I am using Visual Studio 2019 and Parallel Studio 2020.

Thank you very much for your help mecej4. It is much appreciated.

P.S. MKL can be included in Visual Studio as follows: 

Project | Properties | Fortran | Libraries | Use Intel Match kernel Library | Sequential (or Parallel)

module pdata
    implicit none
    integer, parameter :: m=10, n=3 !no. of data pairs, parameters to fit
    double precision, dimension(m) :: xdat = [1,2,3,4,5,6,7,8,9,10]
    double precision, dimension(m) :: ydat = [4,13,28,49,76,109,148,193,244,301]
end module pdata
	
include 'mkl_rci.f90'
program nlfit
    use pdata, only : m, n
    use mkl_rci
    implicit none
    external            func
    integer             rci_request, iter,itrn,st_cr,successful
    integer             iter1, iter2, i, j, info(6)
    double precision    x(n), eps(6), jac_eps, fvec(m), fjac(m, n)
    double precision    rs, r1, r2
    type(handle_tr)     handle
          
    do i = 1, 6
        eps (i) = 1.d-5
    end do
    iter1 = 1000
    iter2 = 100
    rs = 100.d0
    jac_eps = 1.d-8
    x = [0d0,5d0,1.3d0]
    do i = 1, m
        fvec (i) = 0.d0
        do j = 1, n
            fjac (i, j) = 0.d0
        end do
    end do
    if (dtrnlsp_init (handle, n, m, x, eps, iter1, iter2, rs) .ne. tr_success) then
        print *, ' Error in dtrnlsp_init'
        call mkl_free_buffers
        stop 1
    end if
    if (dtrnlsp_check (handle, n, m, fjac, fvec, eps, info) .ne. tr_success) then
        print *, ' Error in dtrnlspbc_init'
        call mkl_free_buffers
        stop 1
    else
        if( info(1) .ne. 0 .or. info(2) .ne. 0 .or. info(3) .ne. 0 .or. info(4) .ne. 0 ) then
            print *, ' Input parameters are not valid'
            call mkl_free_buffers
            stop 1
        end if
    end if
    rci_request = 0
    successful = 0
    itrn = 0
    print *,' itrn  ||f||'
    do while (successful == 0)
        if (dtrnlsp_solve (handle, fvec, fjac, rci_request) .ne. tr_success) then
            print *, ' Error in dtrnlsp_solve'
            call mkl_free_buffers
            stop 1
        end if
        select case (rci_request)
        case (-1, -2, -3, -4, -5, -6)
            successful = 1
        case (1)
            itrn = itrn+1
            call func (m, n, x, fvec)
            print '(i3,es10.2)',itrn,norm2(fvec)
        case (2)
            if (djacobi (func, n, m, fjac, x, jac_eps) .ne. tr_success) then
                print *, ' Error in djacobi'
                call mkl_free_buffers
                stop 1
            end if
        endselect
    end do
    if (dtrnlsp_get (handle, iter, st_cr, r1, r2) .ne. tr_success) then
        print *, ' Error in dtrnlsp_get'
        call mkl_free_buffers
        stop 1
    end if
    if (dtrnlsp_delete (handle) .ne. tr_success) then
        print *, ' Error in dtrnlsp_delete'
        call mkl_free_buffers
        stop 1
    end if

    call mkl_free_buffers
    if (r2 .lt. 1.d-5) then
        PRINT *, 'Best Fit Coefficients'
        print '(3ES13.5)',x(1:n)
    else
        PRINT *, ' Dtrnlsp failed'
	end if
	pause !keeps the console window from disappearing
end program nlfit

subroutine func (m, n, x, f)
    use pdata, only : xdat,ydat
    implicit none
    integer m, n
    double precision x (*), f (*)
    integer i
    
    do i = 1, m
       f(i) = x(1)+x(2)*xdat(i)**x(3)-ydat(i)
    end do
end subroutine func

[Khang_N_Intel]

Since the solution has been given, this thread is now closed.