Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
7002 Discussions

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

Novice
2,666 Views

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.

1 Solution
Novice
2,446 Views

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``````

8 Replies
Moderator
2,603 Views

Hi,

Thanks for reaching out to us.

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

Regards,

Arpita

Honored Contributor III
2,559 Views

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
``````

Novice
2,513 Views

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.

Honored Contributor III
2,507 Views

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.

Novice
2,499 Views

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?

Honored Contributor III
2,490 Views

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.

Novice
2,447 Views

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``````

Employee
2,318 Views

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