使用ctypes在Fortran函数上使用scipy.optimize.minimize的错误结果

Question

我目前有一个Fortran函数,我希望使用SciPy使用Ctypes进行优化.这可能吗？也许我在实施中做错了.例如,假设我有:

cost.f90

module cost_fn
  use iso_c_binding, only: c_float
  implicit none

contains

  function sin_2_cos(x,y) bind(c)
      real(c_float) :: x, y, sin_2_cos
      sin_2_cos = sin(x)**2 * cos(y)
  end function sin_2_cos

end module cost_fn

我编译的:

gfortran -fPIC -shared -g -o cost.so cost.f90

然后尝试找到(本地)最小值:

cost.py

#!/usr/bin/env python

from ctypes import *
import numpy as np
import scipy.optimize as sopt

cost = cdll.LoadLibrary('./cost.so')

cost.sin_2_cos.argtypes = [POINTER(c_float), POINTER(c_float)]
cost.sin_2_cos.restype = c_float

def f2(x):
    return cost.sin_2_cos(c_float(x[0]), c_float(x[1]))
    # return np.sin(x[0])**2 * np.cos(x[1])

# print(f2( [1, 1] ))
# print(f2( [0.5 * np.pi, np.pi] ))

print( sopt.minimize( f2, (1.0, 1.0), options={'disp': True}, tol=1e-8) )

我期望局部最小f2(pi/2,pi)= -1.当我使用cost.sin_2_cos返回值调用f2时,"最小值"仅在(1,1)的初始猜测时给出.如果我使用"Python"返回值调用f2,则optimize会找到正确的最小值.

我已经尝试重新定义sin_2_cos以获取维度(2)数组输入,但是看到了类似的行为.也许我需要用最小化直接调用sin_2_cos(但那么我如何为参数指定c_float)？任何想法都表示赞赏!

编辑:对于下面的评论,请注意两条注释print(f2(...))行产生预期值.因此,我相信通过Python f2函数正确调用Fortran函数.

Answer 1

您的Fortran代码使用单精度浮点值(即32位浮点数).(在C中,float是单精度并且double是双精度.)scipy.optimize.minimize()使用有限差分来近似函数导数的默认方法.也就是说,为了估计导数f'(x0),它计算步长(f(x0+eps) - f(x0))/eps在哪里eps.用于计算有限差分的默认步长大约为1.49e-08.不幸的是,这个值小于值1附近的单精度值的间距.因此,当最小化器eps加1时,结果仍为1.这意味着函数在同一点进行求值,有限差分结果为0这是最小的条件,因此求解器决定它已完成.

求解器选项eps设置有限差分步长.将其设置为大于1.19e-7的值.例如,如果我使用options={'disp': True, 'eps': 2e-6},我会得到一个解决方案.

顺便说一句,您可以使用以下命令找到该值1.19e-7 numpy.finfo():

In [4]: np.finfo(np.float32(1.0)).eps
Out[4]: 1.1920929e-07

如果您method='nelder-mead'在minimize()函数中使用该选项,您还将获得解决方案.该方法不依赖于有限差异.

最后,您可以将Fortran代码转换为使用双精度:

module cost_fn
  use iso_c_binding, only: c_double
  implicit none

contains

  function sin_2_cos(x,y) bind(c)
      real(c_double) :: x, y, sin_2_cos
      sin_2_cos = sin(x)**2 * cos(y)
  end function sin_2_cos

end module cost_fn

然后更改Python代码ctypes.c_double而不是使用ctypes.c_float.