Niepoprawny wynik z scipy.optimize.minimize użyty w funkcji Fortran przy użyciu ctypów

Question

Mam obecnie funkcję Fortran, którą chcę zoptymalizować za pomocą SciPy owijając ją przy pomocy Ctypes. czy to możliwe? Być może popełniłem błąd w mojej implementacji. Na przykład załóżmy, że mam:

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 

że mogę skompilować z:

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

a następnie spróbuj znaleźć (lokalne) minimum z:

kosztów .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)) 

I ex pect lokalne minimum f2 (pi/2, pi) = -1. Kiedy wywołuję f2 z zwracaną wartością cost.sin_2_cos, "minimimum" jest podawane po początkowym odgadnięciu (1, 1). Jeśli wywołasz f2 z wartością zwracaną "Python", optymalizacja znajdzie prawidłowe minimum.

Próbowałem przedefiniować sin_2_cos, aby przyjąć wejście wymiarowe (2), ale widziałem podobne zachowanie. Być może muszę wywołać sin_2_cos bezpośrednio przy zminimalizowaniu (ale w jaki sposób powinienem określić c_float dla argumentów)? Wszelkie przemyślenia są doceniane!

Edytuj: W komentarzu poniżej zauważ, że dwie skomentowane linie print(f2(...)) dają oczekiwane wartości. Dlatego uważam, że funkcja Fortran jest poprawnie wywoływana przez funkcję Python f2.

Answer 1

Twój kod Fortrana wykorzystuje wartości zmiennoprzecinkowe o pojedynczej precyzji (tzn. 32-bitowe zmienne). (W języku C, float jest pojedyncza precyzja, a double jest podwójną precyzją.) Domyślna metoda używana przez scipy.optimize.minimize() używa skończonych różnic do przybliżenia pochodnych funkcji. To znaczy, aby oszacować pochodną f'(x0), oblicza ona (f(x0+eps) - f(x0))/eps, gdzie eps jest wielkością kroku. Domyślny rozmiar kroku, którego używa do obliczenia różnic skończonych, wynosi około 1,49e-08. Niestety, ta wartość jest mniejsza niż odstępy pojedynczych wartości dokładności wokół wartości 1. Tak więc, gdy minimizer dodaje eps do 1, wynik jest nadal 1. Oznacza to, że funkcja jest oceniana w tym samym punkcie, a wynik skończonej różnicy jest 0. Jest to warunek minimum, więc solver decyduje, że to zrobione.

Opcja solver eps ustawia rozmiar kroku różnicowego. Ustaw go na większy niż 1.19e-7. Na przykład, jeśli używam options={'disp': True, 'eps': 2e-6}, otrzymuję rozwiązanie.

Nawiasem mówiąc, można stwierdzić, że wartość, 1.19e-7, za pomocą numpy.finfo():

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

Będziesz również uzyskać rozwiązanie w przypadku korzystania z opcji method='nelder-mead' w funkcji minimize(). Ta metoda nie opiera się na różnicach skończonych.

Wreszcie, można konwertować kod Fortran używać podwójnej precyzji:

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 

Następnie zmień kodu Pythona używać ctypes.c_double zamiast ctypes.c_float.