Discussion :

[romainmill]

Bonjour à tous,

Voici mon problème. J'ai l'équation suivante x=a*cos(b-psi)+c*tan(alpha), avec x et psi des images de 1201x1201.

Du coup, je connais x, psi et alpha et j'aimerai trouvé a, b et c. On m'a dit que je peux utiliser scipy.optimize, mais je vous avouerai que j'ai du mal à m'en servir.

J'ai essayé d'utilisé ce forum mais j'ai du mal à adapter mon problème : http://stackoverflow.com/questions/7...quares-fitting .

Auriez vous des conseils pour m'éclairer, merci d'avance

Merci d'avance pour votre aide

[Julien N]

Salut Romainmill,

Je suis loin d'être un expert en optimisation, mais j'en ai fait un peu. Il existe différentes méthodes plus ou moins complexes pour cela. La plus simple et la plus connue c'est sans doute la méthode des moindres carrés (least mean squares). Scipy vient avec une excellente bibliothèques de fonctions dont leastsq.

La doc est plutôt pas mal, surtout la version PDF. Tu peux notamment y trouver l’exemple sur la fonction de Rosenbrock :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
from scipy.optimize import fmin
 
def rosen(x):
    """The Rosenbrock function"""
    return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
 
x0 = [1.3, 0.7, 0.8, 1.9, 1.2]
xopt = fmin(rosen, x0, xtol=1e-8)
 
print xopt

Ici on cherche a minimiser la somme. Dans ton cas on peut essayer ceci (exemple):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
import numpy as np
from scipy.optimize import fmin
 
def romainmill(x, *args):
    psi = args[0]
    alpha = args[1]
    target = args[2]
    return target - (x[0]*np.cos(x[1]-psi) + x[2]*np.tan(alpha))
 
x0 = [1.0, 2.0, 3.0]
xopt = fmin(romainmill, x0, args=(4.0, 5.0, 0.0))
 
print xopt

La function accepte une liste en entrée (x) qui sont pour toi les trois paramètres à trouver (a, b et c). En plus, une liste d'argument est donnée te permettant de definir alpha, psi et une target (pour minimiser il faut un objectif).

x0 est la liste des valeurs initiales de a, b et c. Les chosir influence le résultat.

C'est vraiment très basique, non linéaire, mais aussi très peu efficace. Tu peux faire bien plus compliqué, mais commence par faire simple.

Ju

[romainmill]

Ok cool merci beaucoup ! c'est plus clair maintenant !

[romainmill]

mais du coup la liste je dois spécifier psi, alpha et target en argument du coup ? et a,b,c aussi ?du genre :

romainmill([a,b,c],[psi,alpha,target]) (cela plante car a, b,c ne sont pas définis

(désolé je commence juste python),

Par ailleurs, quelle est la différence avec la fonction leastsq ?

[Julien N]

Tout à fait, a, b et c ne sont pas définis en tant que tels. Ce sont dans l'ordre les éléments de x, la liste passée en paramètre. Peut-être que le nom x porte à confusion...

En somme on part d'une base, x0. La méthode fmin fonctionne comme une bille sur une surface vallonnée, elle tend à rouler vers le point le bas. Si la surface est très vallonnée, il existe de multiples "vallées" dans lesquelles la bille peut s'arrêter, il existe donc de multiple solution au problème. Je rappel qu'on cherche ici à "minimiser" une quantité, pas à trouver une solution exacte. C'est pourquoi le choix des valeurs initiales est primordiale pour obtenir un résultats intéressant.

Fmin implémente la méthode "simplex":
http://docs.scipy.org/doc/scipy/refe...mize.fmin.html

Leastq fonctionne un peu différemment. Le principe c'est que l'on a une série de valeurs et une loi/fonction que l'on veut coller. En somme trouver une "courbe" (si on est en 2D) passant au plus près des points (valeurs). Pour cela on calcul l'erreur entre notre fonction et les valeurs connues. Cette erreur est élevée au carrée, c'est aussi selon le point de vue une "distance" entre la courbe de la fonction et les points. On somme ces erreurs puis on cherche à minimiser cette somme.

C'est ce que fait excel lorsqu'on cherche à effectuer une régression linéaire.

Bref,

Il ne faut pas faire romainmill([a,b,c],[psi,alpha,target]). Tu peux voir si la fonction romainmill te semble correcte en faisant quelque chose comme ceci:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
a, b, c = 1.0, 2.0, 3.0
psi = 4.0
alpha = 5.0
target = 42.0
args = [psi, alpha, target]
 
err = romainmill([a,b,c], *args)
print err

Fmin va essayer de trouver x[0], x[1], x[2] (qui sont pour toi a, b, c) de telle sorte que x[0]*np.cos(x[1]-psi) + x[2]*np.tan(alpha) soit au plus proche de target = 42.0

Bon tout ça ne résoud pas ton problème. Fmin n'est pas forcément la méthode qui te convient. Est-ce que tu peux nous donner plus de précisions sur tes matrices?

Ouaf, je te laisse digérer ça !

J

[romainmill]

Merci encore pour cette réponse détaillée!

Alors mon problème j'ai deux images de 1201x1201 pixels, identique en apparence mais en fait légérement décalées, je veux évaluer le décalage entre ces deux images.
Celui ci est régis par l'équation que j'ai posté dans le premier message.
En fait a, b et c sont la norme, la direction et la composante verticale du vecteur entre les deux images. Je pensais donc que least square était mieux, car on pouvait évaluer l'erreur entre les deux images (quand l'écart type entre l'image maitre et l'esclave et infèrieure à 2% je dois arrêter les itérations).
Je vais quand même essayé avec fmin pour voir les résultats que cela donne, et après si j'arrive à comprendre leastsq comparer les résultats.

[Julien N]

hmm je dois avouer ne rien y connaitre en traitement d'image. Donc je ne peux pas te dire si tu es sur la bonne voie ou non.


leastsq est plutôt simple à utiliser:
http://docs.scipy.org/doc/scipy/refe...e.leastsq.html

Cette method va grosso modo modifier successivement les valeurs de a, b et c (de façon infinitesimal), et évaluer dans quelles direction "advancer" pour atteindre le minimum. C'est un procéder iteratife. Tu peux grandement améliorer le comportement de la methode en définissant toi meme la matrice jacobienne de ton system.

Tu peux aussi choisir la methode fmin_slsqp, pour le cas où tu souhaiterais contraindre ta recherche d'optimisation. il te faudrait alors définir les trois fonctions suivantes:

1. func(x, *args), c'est la fonction régissant ton décalage
2. eqcons(x, *args), c'est la fonction définissant les égalités (OPTIONEL). Ce sont les contraintes telles que func(x) = target
3. ieqcons(x, *args), c'est la fonction définissant les inégalités (OPTIONEL). Elles sont du type func(x) < val

Il te resterais plus qu'à appeler la methode de scipy:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
res = leastq(func, x0, f_eqcons=eqcons, f_ieqcons=ieqcons,
        bounds=[], iter= 500, iprint=0, full_output=True)

Bon courage,

J

PS: Désolé pour l'orthographe, j'ai un correcteur auto en anglais..

[romainmill]

ok cool.

Je crois que j'ai un peu mieux compris le fonctionnement. Du coup, le probème est que je traite avec des images en effet, et que je ne sais pas comment les mettre en paramètres (x1,x2,x3) renvoie les mm messages que [x1,x2,x3])

Pour fmin j'ai le message d'erreur suivant : (car les trois élèments de args sont des images=matrice de 1201x1201)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-58-70c780962d3d> in <module>()
----> 1 xopt=fmin(coreg.func,x,args)
 
/usr/lib64/python2.7/site-packages/scipy/optimize/optimize.pyc in fmin(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback)
    349             'return_all': retall}
    350 
--> 351     res = _minimize_neldermead(func, x0, args, callback=callback, **opts)
    352     if full_output:
    353         retlist = res['x'], res['fun'], res['nit'], res['nfev'], res['status']
 
/usr/lib64/python2.7/site-packages/scipy/optimize/optimize.pyc in _minimize_neldermead(func, x0, args, callback, xtol, ftol, maxiter, maxfev, disp, return_all, **unknown_options)
    413     if retall:
    414         allvecs = [sim[0]]
--> 415     fsim[0] = func(x0)
    416     nonzdelt = 0.05
    417     zdelt = 0.00025
 
ValueError: setting an array element with a sequence.

et pour la methode avec least square, j'ai :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 
In [60]: scipy.optimize.leastsq(coreg.func,x0,args,full_output=True)---------------------------------------------------------------------------
 
ValueError                                Traceback (most recent call last)
ValueError: object too deep for desired array
---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
<ipython-input-60-a78ee188a7c2> in <module>()
----> 1 scipy.optimize.leastsq(coreg.func,x0,args,full_output=True)
 
/usr/lib64/python2.7/site-packages/scipy/optimize/minpack.pyc in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    362             maxfev = 200*(n + 1)
    363         retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
--> 364                 gtol, maxfev, epsfcn, factor, diag)
    365     else:
    366         if col_deriv:
 
error: Result from function call is not a proper array of floats.

Si je lance un petit whos :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 
Variable      Type        Data/Info
-----------------------------------
alpha         ndarray     1201x1201: 1442401 elems, type `float64`, 11539208 bytes (11 Mb)
args          tuple       n=3
aspect        ndarray     1201x1201: 1442401 elems, type `float64`, 11539208 bytes (11 Mb)
coreg         module      <module 'dem_coreg' from 'dem_coreg.py'>
delta_h       ndarray     1201x1201: 1442401 elems, type `int16`, 2884802 bytes (2 Mb)
denom         ndarray     1201x1201: 1442401 elems, type `float64`, 11539208 bytes (11 Mb)
fmin          function    <function fmin at 0x3e80b18>
fmin_slsqp    function    <function fmin_slsqp at 0x3e8a938>
i             ndarray     1201: 1201 elems, type `float64`, 9608 bytes
leastsq       function    <function leastsq at 0x3e8ad70>
master_dem    ndarray     1201x1201: 1442401 elems, type `int16`, 2884802 bytes (2 Mb)
np            module      <module 'numpy' from '/us<...>ages/numpy/__init__.pyc'>
optimize      module      <module 'scipy.optimize' <...>y/optimize/__init__.pyc'>
pl            module      <module 'pylab' from '/us<...>site-packages/pylab.pyc'>
psi           ndarray     1201x1201: 1442401 elems, type `float64`, 11539208 bytes (11 Mb)
scipy         module      <module 'scipy' from '/us<...>ages/scipy/__init__.pyc'>
x             tuple       n=3
x0            tuple       n=3

[Julien N]

Ok, je pige pas les messages d'erreurs.

Est-ce que tu pourrais poster (si c'est envisageable) un morceau de code me permettant de faire des tests? Dans l'idéal il me faudrait ta fonction coreg.func, la partie faisant appel à scipy, et plus difficile un exemple de matrices. A ce que j'ai compris, tu as deux images rêprésentées par les matrices X et psi. Est-ce que tu crois possible de travailler sur des matrices de plus petites dimensions pour commencer?

Attention à ne pas confondre les variables X (l'une de tes deux images) et x le parametre de tes fonctions. As-tu quelque chose dans ce gout là?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
def func(x, *args):
    psi = args[0]
    alpha = args[1]
    X = args[2]
    return X - (x[0]*np.cos(x[1]-psi) + x[2]*np.tan(alpha))
 
x0 = [1.0, matrice, 3.0]     
args = [psi, alpha, X]   # psi, alpha et X sont des matrices 1201x1201
xopt = fmin(func, x0, *args)
 
print xopt

Attention aux dimensions des élements manipuler. Je suis perdu dans ta formule. X et psi sont toutes deux des matrices 1201x1201 et je suppose que alpha également. Dans ce cas, a et c sont des nombres reels, tandis que b est une matrice (pour pouvoir faire b-psi). Autre hypothèse, a, b et c sont des matrices de meme taille que X et psi, tandis que alpha est un reel...

L'exemple ci-dessus ne marchera pas. La function func() retourne également une matrice de taille 1201x1201. Que veut dire "minimiser" une matrice?

Bref, il faut commencer par éclaircir ce que tu veux minimiser. En 2d, si X, psi et alpha étaient des vecteurs, tu pourrais minimiser la norme de x - (a*cos(b-psi)+c*tan(alpha)). Pour une matrice... Je ne sais pas trop.

Ju

[romainmill]

Ok je vais essayé dêtre plus clair et de te donner les élèments nécessaire pour comprendre.

voici la fonction qui définis mes paramètres que je vais utiliser dans la fonction à minimiser plus tard.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
def coregistration_parameters(image1,image2,band2path,band5path):
    #setting parameters
    master_dem = sf.read_envi_file(image1) #image de 1201X1201
    alpha, psi = slpa.grad2d(master_dem) #calculate slope and aspect of image, alpha et psi deux matrice 1201x1201
    delta_h = pld.plot_difference(image1,image2,-50,50) #différence d'altitude entre les deux images 1201x1201
 
    alpha_barre = np.mean(alpha)
    return delta_h, np.tan(alpha), alpha, alpha_barre, aspect, master_dem

Voici la fonction que je cherche à minimiser : avec x = [a,b,c], a b et c sont des matrices 1201x1201, qui sont les paramètres que je cherche
et args=[aspect,alpha,delta_h]

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
 
 
def func(x, *args):
  psi = args[0]
  alpha = args[1]
  target = args[2]
  return np.array(target - (x[0]*np.cos(x[1]-psi) + x[2]*np.tan(alpha)))

Je t'envoie en mp les trois matrices alpha, aspect et delta_h.

En fait, j'ai compris mes erreurs précédentes en écrivant ce message : j'entrais des valeurs initiales de a,b,c qui était des scalaires. du coup j'ai fait :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
 
In [31]: a0=3.8*np.ones(1442401).reshape((1201,1201))
 
In [33]: b0=1.3*np.ones(1442401).reshape((1201,1201))
 
In [34]: c0=3.1*np.ones(1442401).reshape((1201,1201))

et cette fois c'est le message d'erreur suivant :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
In [39]: scipy.optimize.leastsq(coreg.func, x0, args)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-39-b60529062e3f> in <module>()
----> 1 scipy.optimize.leastsq(coreg.func, x0, args)
 
/usr/lib64/python2.7/site-packages/scipy/optimize/minpack.pyc in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    362             maxfev = 200*(n + 1)
    363         retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
--> 364                 gtol, maxfev, epsfcn, factor, diag)
    365     else:
    366         if col_deriv:
 
ValueError: object too deep for desired array

[Julien N]

Joint plutôt tes matrices sous la forme de fichier texte pour que chacun puisse les exploiter (sauf si ce sont des données sensibles). Pour cela utilise l'icône en forme de trombone.

Pour aisément écrire ces fichiers je te conseils de faire ceci:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

np.savetxt(fname, alpha)

De sorte que je puisse recréer la matrice en utilisant la methode inverse:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

np.loadtxt(fname)

Ju

[romainmill]

impossible via ce site depuis mon ordinateur. L'interface gérer les pièces jointes n'ajoute pas les fichiers txt que je demande. Par mail si vous préférez je peux aussi pour ceux qui veulent tester le code

[romainmill]

a,b et c sont bien des réels en fait. Il est possible de faire un nombre réel-une matrice. je fixe (a=3,b=1,c=0)

[Julien N]

Demain je suis en déplacement, je pourrais pas trop regarder. Peut-être que d'ici là quelqu'un t'aura apporté une réponse.

Si tu peux faire a- X, c'est donc que tu fais a-Xij pour chaque i, j de ta matrice. Pourquoi pas.

[romainmill]

Ok pas de soucis ! je vais continuer à plancher dessus demain aussi

[romainmill]

Pour essayer d'éclaircir encore une fois les choses.

J'essaye en fait de trouver les paramètres qui vont fitter au mieux mes jeux de données, (voir figure jointe).


La fonction qui régis ce jeu de donnée est x[0]*np.cos(x[1]-psi) + x[2]). Je cherche donc à minimiser : target-(x[0]*np.cos(x[1]-psi) + x[2]) pour trouver les meilleurs x[1],x[0],x[2].

Mes données sont alors psi et target que je range dans args
Mes paramètres initiaux sont rangés dans x, et sont a0,b0,c0

je lance donc la commande :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
leastsq(coreg.func,x,args)

Encore une fois, je veux tirés de cette minimisation mes trois paramètre a,b,c qui ajuste au mieux mes données.


Pour avoir une idée des objets que je manipule :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
 
args           tuple          n=2
a0             float          3.8
b0             float          1.3
c0             int            0
target         ndarray        1201x1201: 1442401 elems, type `float64`, 11539208 bytes (11 Mb)
x              ndarray        3: 3 elems, type `float64`, 24 bytes
psi            ndarray        1201x1201: 1442401 elems, type `float64`, 11539208 bytes (11 Mb)

voici le message d'erreur reçu, et donc le soucis que je n'arrive pas à élucider...

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 
In [104]: leastsq(coreg.func,x_v,args_v)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
ValueError: object too deep for desired array
---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
<ipython-input-104-7a02a43bc618> in <module>()
----> 1 leastsq(coreg.func,x_v,args_v)
 
/usr/lib64/python2.7/site-packages/scipy/optimize/minpack.pyc in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    362             maxfev = 200*(n + 1)
    363         retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
--> 364                 gtol, maxfev, epsfcn, factor, diag)
    365     else:
    366         if col_deriv:
 
error: Result from function call is not a proper array of floats.

[romainmill]

Bon, je ne sais pas si j'ai réussi à régler le problème, mais il n'affiche plus le message d'erreur. En fait ma fonction qui au départ était définie tel quel :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
def func(x, *args):
  psi = args[0]
  target = args[1]
  return np.array(target - (x[0]*np.cos(x[1]-psi) + x[2]))

s'est transformée en :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
def func(x, *args):
  psi = args[0].ravel()
  target = args[1].ravel()
  return np.array(target - (x[0]*np.cos(x[1]-psi) + x[2]))

et le message d'erreur ne s'affiche plus. Par contre il me renvoie la chose suivante :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
In [96]: leastsq(coreg.func,x0,args,full_output=1)
Out[96]: 
(array([ 3.,  1.,  0.]),
 None,
 {'fjac': array([[ nan,  nan,  nan, ...,  nan,  nan,  nan],
       [ nan,  nan,  nan, ...,  nan,  nan,  nan],
       [ nan,  nan,  nan, ...,  nan,  nan,  nan]]),
  'fvec': array([-2.3222515 , -2.02262394, -1.86777113, ...,  2.02609317,
        2.35137464,  1.63630145]),
  'ipvt': array([1, 2, 3], dtype=int32),
  'nfev': 4,
  'qtf': array([ nan,  nan,  nan])},
 'The cosine of the angle between func(x) and any column of the\n  Jacobian is at most 0.000000 in absolute value',
 4)

C'est à dire exactement les mêmes paramètre que ce que j'ai en entrée (et ça ne fit pas terrible)

[Julien N]

Salut,

Je n'ai pas eu le temps de m'y mettre, désolé.

Une chose, pourquoi transformer tes matrices 2D en vecteur 1D? Est-ce que ce faisant ta function d'optimisation reste correcte? Est-ce toujours cela que tu veux minimizer?

Le problem visiblement c'est la matrice jacobienne ne peut être calculée. Essaye ceci:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
def func(x, *args):
  psi = args[0].ravel()
  target = args[1].ravel()
  return np.linalg.norm(target - (x[0]*np.cos(x[1]-psi) + x[2]))

J

[romainmill]

Salut ! =)

Alors sans transformer ma matrice 2D en 1D, la fonction d'optimisation affiche le message précédent (not a proper array of float). Après je peux toujours en sortie faire un reshape pour retransformer en matrice ce qui sort d ema fonction d'optimisation?

Sinon j'ai essayé ta suggestion et voilà ce qu'il se passe :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
In [13]: leastsq(coreg.func,x0,args,full_output=1)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-13-8e7be47ef916> in <module>()
----> 1 leastsq(coreg.func,x0,args,full_output=1)
 
/usr/lib64/python2.7/site-packages/scipy/optimize/minpack.pyc in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    355     m = shape[0]
    356     if n > m:
--> 357         raise TypeError('Improper input: N=%s must not exceed M=%s' % (n,m))
    358     if epsfcn is None:
    359         epsfcn = sqrt(finfo(dtype).eps)
 
TypeError: Improper input: N=3 must not exceed M=1

ce qui est surement du au fait que (si j'utilise mes valeurs initial pour x):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
 
In [24]: np.linalg.norm(target - (x[0]*np.cos(x[1]-psi) + x[2]))
Out[24]: nan

[Julien N]

Il y a quelque chose qui ne va pas.

Voici ce que j'ai chez moi:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# -*- coding: utf-8 -*-
 
import numpy as np
from scipy.optimize import fmin_slsqp
 
# Load matrices
alpha = np.loadtxt('alpha.txt')
aspect = np.loadtxt('aspect.txt')
delta_h = np.loadtxt('delta_h.txt')
 
 
def func(x, *args):
    psi = args[0].ravel()
    target = args[1].ravel()
    return np.linalg.norm(target - (x[0]*np.cos(x[1]-psi) + x[2]))
 
 
if __name__ == '__main__':
 
    x0 = [-8.0, 8.0, -3.0]
    args = [aspect, delta_h]
 
    # Print start value
    print func(x0, *args)
 
    # Least Square Method
    xopt = fmin_slsqp(func, x0, args=args, iter= 500, iprint=0, full_output=True) 
    print xopt

La version de Scipy que j'utilise est la 10.1 (pas le choix, c'est celle qui est fournie au bureau). La méthode fmin_slsqp a changée de nom depuis pour leastsq .

Je me plante peut-être dans l'utilisation de tes matrices, je suis perdu entre les différentes appellations: target, psi, aspect, delta_h... Et tu n'utilise plus non plus la matrice alpha dans tes derniers posts.

Bref, première chose, appeler la fonction avec les valeurs de départs fonctionne:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
x0 = [-8.0, 8.0, -3.0]
args = [aspect, delta_h]
print func(x0, *args)
>>> 19506.37

Mieux, la méthode d'optimisation tourne et sort quelque chose:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
xopt = fmin_slsqp(func, x0, args=args, iter= 500, iprint=0, full_output=True) 
print xopt
>>> [[-7.6110908956373535, 8.0177891995050725, -2.8127125487092268], 19501.816783155
31, 8, 0, 'Optimization terminated successfully.']

Je ne crois pas aux coincidences. Si les valeurs trouvées sont si proches de mes valeurs de depart, c'est qu'il y a de très fortes chances qu'il existe de multiple minimums. C'est pas gagné.

Une dernière chose, ici on fait une norme, c'est pour avoir une valeur numérique à minimiser. J'ai du mal à voir comment ordonner des vecteurs, et encore moins des matrices. La function d'optimisation actuelle est donc peut-être mauvaise, mais ça c'est à toi de voir.

Ju