使用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 

，然後嘗試找到一個（本地）最低配：

成本.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)) 

我ex將局部最小值f2（pi/2，pi）= -1。當我用cost.sin_2_cos返回值調用f2時，「minimimum」僅在（1,1）的初始猜測中給出。如果我用「Python」返回值調用f2，則優化會找到正確的最小值。

我試過重新定義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 

您還可以，如果你在minimize()功能使用選項method='nelder-mead'得到解決。該方法不依賴於有限差異。

最後，你可以轉換的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。