<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 2A.data - Calcul Matriciel, Optimisation
numpy arrays sont la première chose à considérer pour accélérer un algorithme. Les matrices sont présentes dans la plupart des algorithmes et numpy optimise les opérations qui s'y rapporte. Ce notebook est un parcours en diagonal.
Step1:  Numpy arrays
La convention d'import classique de numpy est la suivante
Step2:  Creation d'un array
Step3:  On peut indiquer explicitement le dtype lors de la création de l'array. Sinon, Numpy sélectionne automatiquement le dtype.
Numpy ajoute un grand nombre de dtype à ceux de Python. Allez jeter un oeil à la liste.
Step4:  Assigner un float dans un array de type int va caster le float en int, et ne modifie pas le dtype de l'array.
Step5:  On peut forcer le casting dans un autre type avec astype 
Step6:  A partir d'une liste de listes, on obtient un array bi-dimmensionnel.
On peut le transposer ou encore l'aplatir en un array 1d
Step7:  Indexation, Slicing, Fancy indexing
Step8:  L'indexation des array multidimensionnels fonctionne avec des tuples.
La syntaxe '
Step9:  Si on utilise pas un couple sur un array 2d on récupère un array 1d
Step10:  On peut aussi utiliser l'indexation par un array (ou une liste python) de booléens ou d'entiers (un mask). Cela s'appelle le fancy indexing. Un mask d'entiers permet de désigner les éléments que l'on souhaite extraire via la liste de leurs indices, on peut aussi répéter l'indice d'un élément pour répéter l'élement dans l'array que l'on extrait.
Step11:  Pourquoi parle-t-on de fancy indexing? Essayez d'indexer des listes python de la même manière...
Step12:  View contre Copy
Créons un array $d$. En plus de renvoyer directement un array, la fonction arange permet aussi d'utiliser un step flottant. (Essayer avec le range de python pour voir)
Step13:  Un point important est que l'on ne recopie pas un array lorsqu'on effectue une assignation ou un slicing d'un array.
On travaille dans ce cas avec une View sur l'array d'origine (shallow copy). Toute modification sur la View affecte l'array d'origine.
Dans l'exemple qui suit, $e$ est une view sur $d$. Lorsqu'on modifie $e$, $d$ aussi est modifié. (Remarquez au passage que numpy fournit quelques constantes bien pratiques....)
Step14:  Si on ne veut pas modifier $d$ indirectement, il faut travailler sur une copie de $d$ (deep copy).
Step15:  Ce point est important car source classique d'erreurs silencieuses
Step16:  Manipulation de shape
La méthode reshape permet de changer la forme de l'array. Il existe de nombreuses manipulations possibles.
On précise à reshape la forme souhaitée
Step17:  Par défaut, reshape utilise l'énumération dans l'ordre du langage C (aussi appelé "row first" ), on peut préciser que l'on souhaite utiliser l'ordre de Fortran ("column first"). Ceux qui connaissent Matlab et R sont habitués à l'ordre "column-first". Voir l'article wikipedia
Step18:  On peut utiliser -1 sur une dimension, cela sert de joker
Step19:  On peut aussi concatener ou stacker horizontalement/verticalement différents arrays.
Step20:  Exercice 1
Step21:  Manipulation et Opérations sur les arrays
Il existe un très grand nombre de routines pour manipuler les arrays numpy
Step22:  Les opérations arithmétiques avec les scalaires, ou entre arrays s'effectuent élément par élément.
Lorsque le dtype n'est pas le même ($a$ contient des float, $b$ contient des int), numpy adopte le type le plus "grand" (au sens de l'inclusion).
Step23:  Remarquez que la division par zéro ne provoque pas d'erreur mais introduit la valeur inf 
Step24:  Broadcasting
Que se passe-t-il si les dimensions sont différentes?
Step25:  L'opération précédente fonctionne car numpy effectue ce qu'on appelle un broadcasting de c 
Step26:  Par contre, il peut parfois être utile de préciser la dimension sur laquelle on souhaite broadcaster, on ajoute alors explicitement une dimension 
Step27:  Réductions
On parle de réductions lorsque l'opération réduit la dimension de l'array.
Il en existe un grand nombre. Elles existent souvent sous forme de fonction de numpy ou de méthodes d'un array numpy.
On n'en présente que quelques unes, mais le principe est le même 
Step28:  Algèbre linéaire
Vous avez un éventail de fonctions pour faire de l'algèbre linéaire dans numpy ou dans scipy.
Cela peut vous servir si vous cherchez à faire une décomposition matricielle particulière (LU, QR, SVD,...), si vous vous intéressez aux valeurs propres d'une matrice, etc...
Exemples simples
Commençons par construire deux arrays 2d correspondant à une matrice triangulaire inférieure et une matrice diagonale 
Step29:  On a vu que les multiplications entre array s'effectuaient élément par élement.
Si l'on souhaite faire des multiplications matricielles, il faut utiliser la fonction dot. La version 3.5 introduit un nouvel opérateur @ qui désigne explicitement la multiplication matricielle.
Step30:  On peut calculer l'inverse ou le déterminant de $A$
Step31:  ... résoudre des systèmes d'equations linéaires du type $Ax = b$...
Step32:  ... ou encore obtenir les valeurs propres de $A$.
Step33:  Numpy Matrix
Matrix est une sous classe spécialisée pour le calcul matriciel. Il s'agit d'un array numpy 2d qui conserve sa dimension 2d à travers les opérations. Pensez aux différences que cela implique...
On peut les construire classiquement depuis les array ou les objets pythons, ou via une string à la Matlab ( où les points virgules indiquent les lignes).
Step34:  Matrix surcharge par ailleurs les opérateurs * et ** pour remplacer les opérations élément par élément par les opérations matricielles.
Enfin, une Matrix possède des attributs supplémentaires. Notamment, Matrix.I qui désigne l'inverse, Matrix.A l'array de base. 
Il est probable que cela évolue puisque Python 3.5 a introduit le symbol @ pour la multiplication matricielle.
Step35:  La syntaxe est plus légère pour effectuer du calcul matriciel
Step36:  Génération de nombres aléatoires et statistiques
Le module numpy.random apporte à python la possibilité de générer un échantillon de taille $n$ directement, alors que le module natif de python ne produit des tirages que un par un. Le module numpy.random est donc bien plus efficace si on veut tirer des échantillon conséquents. Par ailleurs, scipy.stats fournit des méthodes pour un très grand nombre de distributions et quelques fonctions classiques de statistiques. 
Par exemple, on peut obtenir un array 4x3 de tirages gaussiens standard (soit en utilisant randn ou normal)
Step37:  Pour se convaincre que numpy.random est plus efficace que le module random de base de python. On effectue un grand nombre de tirages gaussiens standard, en python pur et via numpy.
Step38:  Exercice 2 
Step39:  On va ensuite chercher un exemple dans la gallerie matplotlib pour la représenter
Step40:  On voit que le minimum se trouve près de $[-\frac{1}{2}, 0]$. On va utiliser ce point pour initialiser l'optimisation.
On va tester différentes méthodes et comparer les sorties obtenues.
Step41:  On trouve un minimum à $-0.4052$ en $[-0.669, 0.000]$ pour toutes les méthodes qui convergent. Notez le message de sortie de 'CG' qui signifie que le gradient ne varie plus assez. Personnellement, je ne trouve pas ce message de sortie très clair. Le point trouvé est bien l'optimum cherché pourtant. Notez aussi le nombre d'évaluations de la fonction (nfev) pour chaque méthode, et le nombre d'évaluation de gradient (njev) pour les méthodes qui reposent sur un calcul de gradient.
Remarquez aussi que si on relance Anneal plusieurs fois, on n'est pas assuré d'obtenir la même solution, puisqu'il s'agit d'une métaheuristique.
Step42:  On va évaluer le temps de calcul nécessaire à chaque méthode.
Step43:  On peut aussi fournir des arguments supplémentaires à la fonction qu'on optimise. Par exemple, les données lorsque vous maximisez une log-vraissemblance. En voici un exemple
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from jyquickhelper import add_notebook_menu
add_notebook_menu()
Explanation: 2A.data - Calcul Matriciel, Optimisation
numpy arrays sont la première chose à considérer pour accélérer un algorithme. Les matrices sont présentes dans la plupart des algorithmes et numpy optimise les opérations qui s'y rapporte. Ce notebook est un parcours en diagonal.
End of explanation
import numpy as np
Explanation: Numpy arrays
La convention d'import classique de numpy est la suivante:
End of explanation
l = [1, 42, 18 ]
a = np.array(l)
print(a)
print(a.dtype)
print(a.ndim)
print(a.shape)
print(a.size)
a
Explanation: Creation d'un array: notion de datatype, et dimensions
On part d'une liste python contenant des entiers. On peut créer un array numpy à partir de cette liste. 
Cet array possède des attributs indiquant le data type, le nombre de dimensions de l'array, etc...
End of explanation
b = np.array(l, dtype=float)
print(b)
print(b.dtype)
l[0] = 1.0
bb = np.array(l)
print(bb)
print(bb.dtype)
Explanation: On peut indiquer explicitement le dtype lors de la création de l'array. Sinon, Numpy sélectionne automatiquement le dtype.
Numpy ajoute un grand nombre de dtype à ceux de Python. Allez jeter un oeil à la liste.
End of explanation
a[0] = 2.5 
a
Explanation: Assigner un float dans un array de type int va caster le float en int, et ne modifie pas le dtype de l'array.
End of explanation
aa = a.astype(float)
aa[0] = 2.5
aa
Explanation: On peut forcer le casting dans un autre type avec astype :
End of explanation
c = np.array([range(5), range(5,10), range(5)])
print(c)
print("ndim:{}".format(c.ndim))
print("shape:{}".format(c.shape))
print(c.transpose()) #same as c.T
print("shape transposed:{}".format(c.T.shape))
print(c.flatten())
print("ndim flattened:{}".format(c.flatten().ndim))
Explanation: A partir d'une liste de listes, on obtient un array bi-dimmensionnel.
On peut le transposer ou encore l'aplatir en un array 1d
End of explanation
print(c)
Explanation: Indexation, Slicing, Fancy indexing
End of explanation
print(c[1,3])
print(c[1,:3])
print(c[:,4])
Explanation: L'indexation des array multidimensionnels fonctionne avec des tuples.
La syntaxe ':' permet d'obtenir tous les éléments de la dimension.
End of explanation
print(c[1], c[1].shape)
print(c[1][:3])
Explanation: Si on utilise pas un couple sur un array 2d on récupère un array 1d
End of explanation
ar = np.arange(1,10) #arange est l'equivalent de range mais retourne un numpy array
print('ar = ',ar)
idx = np.array([1, 4, 3, 2, 1, 7, 3])
print('idx = ',idx)
print("ar[idx] =", ar[idx])
print('######')
idx_bool = np.ones(ar.shape, dtype=bool)
idx_bool[idx] = False
print('idx_bool = ', idx_bool)
print('ar[idx_bool] = ', ar[idx_bool])
print('######', 'Que se passe-t-il dans chacun des cas suivants?', '######' )
try:
    print('ar[np.array([True, True, False, True])] = ', ar[np.array([True, True, False, True])])
except Exception as e:
    # l'expression ar[[True, True, False, True]] déclenche une erreur depuis numpy 1.13
    print("Erreur", e)
Explanation: On peut aussi utiliser l'indexation par un array (ou une liste python) de booléens ou d'entiers (un mask). Cela s'appelle le fancy indexing. Un mask d'entiers permet de désigner les éléments que l'on souhaite extraire via la liste de leurs indices, on peut aussi répéter l'indice d'un élément pour répéter l'élement dans l'array que l'on extrait.
End of explanation
list_python = range(10)
list_python[[True, True, False, True]] # déclenche une exception
list_python[[2, 3, 2, 7]]  # déclenche une exception
Explanation: Pourquoi parle-t-on de fancy indexing? Essayez d'indexer des listes python de la même manière...
End of explanation
d = np.arange(1, 6, 0.5)
d
Explanation: View contre Copy
Créons un array $d$. En plus de renvoyer directement un array, la fonction arange permet aussi d'utiliser un step flottant. (Essayer avec le range de python pour voir)
End of explanation
e = d
e[[0,2, 4]] = - np.pi
e
d
Explanation: Un point important est que l'on ne recopie pas un array lorsqu'on effectue une assignation ou un slicing d'un array.
On travaille dans ce cas avec une View sur l'array d'origine (shallow copy). Toute modification sur la View affecte l'array d'origine.
Dans l'exemple qui suit, $e$ est une view sur $d$. Lorsqu'on modifie $e$, $d$ aussi est modifié. (Remarquez au passage que numpy fournit quelques constantes bien pratiques....)
End of explanation
d = np.linspace(1,5.5,10) #Question subsidiaire: en quoi est-ce différent de np.arange avec un step float?
f = d.copy()
f[:4] = -np.e #il s'agit du nombre d'euler, pas de l'array e ;)
print(f)
print(d)
Explanation: Si on ne veut pas modifier $d$ indirectement, il faut travailler sur une copie de $d$ (deep copy).
End of explanation
print('d = ',d)
slice_of_d = d[2:5]
print('\nslice_of_d = ', slice_of_d)
slice_of_d[0] = np.nan
print('\nd = ', d)
mask = np.array([2, 3, 4])
fancy_indexed_subarray = d[mask]
print('\nfancy_indexed_subarray = ', fancy_indexed_subarray)
fancy_indexed_subarray[0] = -2
print('\nd = ', d)
Explanation: Ce point est important car source classique d'erreurs silencieuses: les erreurs les plus vicieuses car l'output sera faux mais python ne râlera pas...
Il faut un peu de temps pour s'habituer mais on finit par savoir de manière naturelle quand on travaille sur une view, quand on a besoin de faire une copie explicitement, etc... En tout cas, vérifiez vos sorties, faites des tests de cohérence, cela ne nuit jamais.
Retenez par exemple que le slicing vous renvoie une view sur l'array, alors que le fancy indexing effectue une copie.
(Au passage, remarquez le NaN (=NotaNumber) déjà introduit lors de la séance 1 sur pandas qui est un module basé sur numpy)
End of explanation
g = np.arange(12)
print(g)
g.reshape((4,3))
Explanation: Manipulation de shape
La méthode reshape permet de changer la forme de l'array. Il existe de nombreuses manipulations possibles.
On précise à reshape la forme souhaitée: par un entier si on veut un array 1d de cette longueur, ou un couple pour un array 2d de cette forme.
End of explanation
g.reshape((4,3), order='F')
Explanation: Par défaut, reshape utilise l'énumération dans l'ordre du langage C (aussi appelé "row first" ), on peut préciser que l'on souhaite utiliser l'ordre de Fortran ("column first"). Ceux qui connaissent Matlab et R sont habitués à l'ordre "column-first". Voir l'article wikipedia
End of explanation
np.zeros_like(g)
np.ones_like(g)
Explanation: On peut utiliser -1 sur une dimension, cela sert de joker: numpy infère la dimension nécessaire ! On peut créer directement des matrices de 0 et de 1 à la dimension d'un autre array.
End of explanation
np.concatenate((g, np.zeros_like(g))) #Attention à la syntaxe: le type d'entrée est un tuple!
gmat = g.reshape((1, len(g)))
np.concatenate((gmat, np.ones_like(gmat)), axis=0)
np.concatenate((gmat, np.ones_like(gmat)), axis=1)
np.hstack((g, g))
np.vstack((g,g))
Explanation: On peut aussi concatener ou stacker horizontalement/verticalement différents arrays.
End of explanation
#Exo1a-1:
#Exo1a-2:
#Exo1B:
#Exo1C:
Explanation: Exercice 1: Echiquier et Crible d'Erathosthène
Exercice 1-A Echiquier: Créer une matrice échiquier (des 1 et des 0 alternés) de taille 8x8, de deux façons différentes
en vous servant de slices 
en vous servant de la fonction tile
Exercice 1-B Piège lors d'une extraction 2d:
Définir la matrice $M = \left(\begin{array}{ccccc} 1 & 5 & 9 & 13 & 17 \ 2 & 6 & 10 & 14 & 18 \ 3 & 7 & 11 & 15 & 19 \ 4 & 8 & 12 & 16 & 20 \ \end{array}\right)$
En extraire la matrice $\left(\begin{array}{ccc} 6 & 18 & 10 \ 7 & 19 & 11 \ 5 & 17 & 9 \ \end{array}\right)$
Exercice 1-C Crible d'Erathosthène: On souhaite implémenter un crible d'Erathosthène pour trouver les nombres premiers inférieurs à $N=1000$.
partir d'un array de booléens de taille N+1, tous égaux à True.
Mettre 0 et 1 à False car ils ne sont pas premiers
pour chaque entier $k$ entre 2 et $\sqrt{N}$: 
si $k$ est premier: on passe ses multiples (entre $k^2$ et $N$) à False
on print la liste des entiers premiers
End of explanation
a = np.ones((3,2))
b = np.arange(6).reshape(a.shape)
print(a)
b
Explanation: Manipulation et Opérations sur les arrays
Il existe un très grand nombre de routines pour manipuler les arrays numpy: 
Vous trouverez sans doute utiles les pages spécifiques aux routines de stats ou de maths
Opérations élément par élément
On déclare $a$ et $b$ sur lesquelles nous allons illustrer quelques opérations
End of explanation
print( (a + b)**2 )
print( np.abs( 3*a - b ) )
f = lambda x: np.exp(x-1)
print( f(b) )
Explanation: Les opérations arithmétiques avec les scalaires, ou entre arrays s'effectuent élément par élément.
Lorsque le dtype n'est pas le même ($a$ contient des float, $b$ contient des int), numpy adopte le type le plus "grand" (au sens de l'inclusion).
End of explanation
b
1/b
Explanation: Remarquez que la division par zéro ne provoque pas d'erreur mais introduit la valeur inf :
End of explanation
c = np.ones(6)
c
b+c   # déclenche une exception
c = np.arange(3).reshape((3,1))
print(b,c, sep='\n')
b+c
Explanation: Broadcasting
Que se passe-t-il si les dimensions sont différentes?
End of explanation
a = np.zeros((3,3))
a[:,0] = -1
b = np.array(range(3))
print(a + b)
Explanation: L'opération précédente fonctionne car numpy effectue ce qu'on appelle un broadcasting de c : une dimension étant commune, tout se passe comme si on dupliquait c sur la dimension non-partagée avec b. Vous trouverez une explication visuelle simple ici :
End of explanation
print(b.shape)
print(b[:,np.newaxis].shape) 
print(b[np.newaxis,:].shape)
print( a + b[np.newaxis,:] )
print( a + b[:,np.newaxis] )
print(b[:,np.newaxis]+b[np.newaxis,:])
print(b + b)
Explanation: Par contre, il peut parfois être utile de préciser la dimension sur laquelle on souhaite broadcaster, on ajoute alors explicitement une dimension :
End of explanation
c = np.arange(10).reshape((2,-1)) #Note: -1 is a joker! 
print(c)
print(c.sum())
print(c.sum(axis=0))
print(np.sum(c, axis=1))
print(np.all(c[0] < c[1]))
print(c.min(), c.max())
print(c.min(axis=1))
Explanation: Réductions
On parle de réductions lorsque l'opération réduit la dimension de l'array.
Il en existe un grand nombre. Elles existent souvent sous forme de fonction de numpy ou de méthodes d'un array numpy.
On n'en présente que quelques unes, mais le principe est le même : par défaut elles opèrent sur toutes les dimensions, mais on peut via l'argument axis préciser la dimension selon laquelle on souhaite effectuer la réduction.
End of explanation
A = np.tril(np.ones((3,3)))
A
b = np.diag([1,2, 3])
b
Explanation: Algèbre linéaire
Vous avez un éventail de fonctions pour faire de l'algèbre linéaire dans numpy ou dans scipy.
Cela peut vous servir si vous cherchez à faire une décomposition matricielle particulière (LU, QR, SVD,...), si vous vous intéressez aux valeurs propres d'une matrice, etc...
Exemples simples
Commençons par construire deux arrays 2d correspondant à une matrice triangulaire inférieure et une matrice diagonale :
End of explanation
print(A.dot(b))
print(A*b)
print(A.dot(A))
Explanation: On a vu que les multiplications entre array s'effectuaient élément par élement.
Si l'on souhaite faire des multiplications matricielles, il faut utiliser la fonction dot. La version 3.5 introduit un nouvel opérateur @ qui désigne explicitement la multiplication matricielle.
End of explanation
print(np.linalg.det(A))
inv_A = np.linalg.inv(A)
print(inv_A)
print(inv_A.dot(A))
Explanation: On peut calculer l'inverse ou le déterminant de $A$
End of explanation
x = np.linalg.solve(A, np.diag(b))
print(np.diag(b))
print(x)
print(A.dot(x))
Explanation: ... résoudre des systèmes d'equations linéaires du type $Ax = b$...
End of explanation
np.linalg.eig(A)
np.linalg.eigvals(A)
Explanation: ... ou encore obtenir les valeurs propres de $A$.
End of explanation
m = np.matrix(' 1 2 3; 4 5 6; 7 8 9')
a = np.arange(1,10).reshape((3,3))
print(m)
print(a)
print(m[0], a[0])
print(m[0].shape, a[0].shape)
Explanation: Numpy Matrix
Matrix est une sous classe spécialisée pour le calcul matriciel. Il s'agit d'un array numpy 2d qui conserve sa dimension 2d à travers les opérations. Pensez aux différences que cela implique...
On peut les construire classiquement depuis les array ou les objets pythons, ou via une string à la Matlab ( où les points virgules indiquent les lignes).
End of explanation
m * m
a * a
m * a # La priorité des matrix est plus importantes que celles des arrays
print(m**2)
print(a**2)
Explanation: Matrix surcharge par ailleurs les opérateurs * et ** pour remplacer les opérations élément par élément par les opérations matricielles.
Enfin, une Matrix possède des attributs supplémentaires. Notamment, Matrix.I qui désigne l'inverse, Matrix.A l'array de base. 
Il est probable que cela évolue puisque Python 3.5 a introduit le symbol @ pour la multiplication matricielle.
End of explanation
m[0,0]= -1
print("det", np.linalg.det(m), "rank",np.linalg.matrix_rank(m))
print(m.I*m)
a[0,0] = -1
print("det", np.linalg.det(a), "rank",np.linalg.matrix_rank(a))
print(a.dot(np.linalg.inv(a)))
Explanation: La syntaxe est plus légère pour effectuer du calcul matriciel
End of explanation
np.random.randn(4,3)
Explanation: Génération de nombres aléatoires et statistiques
Le module numpy.random apporte à python la possibilité de générer un échantillon de taille $n$ directement, alors que le module natif de python ne produit des tirages que un par un. Le module numpy.random est donc bien plus efficace si on veut tirer des échantillon conséquents. Par ailleurs, scipy.stats fournit des méthodes pour un très grand nombre de distributions et quelques fonctions classiques de statistiques. 
Par exemple, on peut obtenir un array 4x3 de tirages gaussiens standard (soit en utilisant randn ou normal):
End of explanation
N = int(1e7)
from random import normalvariate
%timeit [normalvariate(0,1) for _ in range(N)]
%timeit np.random.randn(N)
Explanation: Pour se convaincre que numpy.random est plus efficace que le module random de base de python. On effectue un grand nombre de tirages gaussiens standard, en python pur et via numpy.
End of explanation
def bowl_peak(x,y): 
    return x*np.exp(-x**2-y**2)+(x**2+y**2)/20
Explanation: Exercice 2 : marches aléatoires
Simulez (en une seule fois!) 10000 marches aléatoires de taille 1000, partant de 0 et de pas +1 ou -1 équiprobables 
Faites un graphe représentant la racine de la moyenne des carrés des positions (=cumul des pas à un instant donné) en fonction du temps
Quels sont les amplitudes maximales et minimales atteintes parmi l'ensemble des marches aléatoires?
Combien de marches s'éloigne de plus de 50 de l'origine?
Parmi celles qui le font, quelle est la moyenne des temps de passage (i.e. le premier moment où ces marches dépassent +/-50)?
Vous aurez peut-être besoin des fonctions suivantes: np.abs, np.mean, np.max, np.where, np.argmax, np.any, np.cumsum, np.random.randint.
Exercice 3 : retrouver la série aléatoire à partir des marches aléatoires
L'exercice précédent montre comment générer une marche aléatoire à partir d'une série temporelle aléatoire. Comment retrouver la série initiale à partir de la marche aléatoire ?
Optimisation avec scipy
Le module scipy.optimize fournit un panel de méthodes d'optimisation. En fonction du problème que vous souhaitez résoudre, il vous faut choisir la méthode adéquate. Je vous conseille vivement la lecture de ce tutoriel sur l'optimisation numérique, écrit par Gaël Varoquaux. 
Récemment, l'ensemble des solvers ont été regroupés sous deux interfaces, même si on peut toujours faire appel à chaque solver directement, ce qui n'est pas conseillé car les entrées sorties ne sont pas normalisées (par contre vous devrez sans doute aller voir l'aide de chaque méthode pour vous en servir): 
Pour minimiser une fonction scalaire d'une ou plusieurs variables:scipy.optimize.minimize
Pour minimiser une fonction scalaire d'une variable uniquement:scipy.optimize.minimize_scalar
Vous obtiendrez en sortie un objet de type scipy.optimize.OptimizeResult.
Dans la suite, je développe un petit exemple inspiré du tutoriel de la toolbox d'optimisation de Matlab. Par ailleurs, la documentation de cette toolbox est plutôt claire et peut toujours vous servir lorsque que vous avez besoin de vous rafraichir la mémoire sur l'optimisation numérique.
On commence par définir la fonction bowl_peak
End of explanation
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm #colormaps
min_val = -2
max_val = 2
fig = plt.figure()
ax = fig.gca(projection='3d')
x_axis = np.linspace(min_val,max_val,100)
y_axis = np.linspace(min_val,max_val,100)
X, Y = np.meshgrid(x_axis, y_axis, copy=False, indexing='xy')
Z = bowl_peak(X,Y)
#X, Y, Z = axes3d.get_test_data(0.05)
ax.plot_surface(X, Y, Z, rstride=5, cstride=5, alpha=0.2)
cset = ax.contour(X, Y, Z, zdir='z', offset=-0.5, cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='x', offset=min_val, cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='y', offset=max_val, cmap=cm.coolwarm)
ax.set_xlabel('X')
ax.set_xlim(min_val, max_val)
ax.set_ylabel('Y')
ax.set_ylim(min_val, max_val)
ax.set_zlabel('Z')
ax.set_zlim(-0.5, 0.5)
Explanation: On va ensuite chercher un exemple dans la gallerie matplotlib pour la représenter: contour3d_demo3. On modifie légèrement le code pour l'utiliser avec bowl_peak
End of explanation
from scipy import optimize
x0 = np.array([-0.5, 0])
fun = lambda x: bowl_peak(x[0],x[1])
methods = [ 'Nelder-Mead', 'CG', 'BFGS', 'Powell', 'COBYLA', 'L-BFGS-B' ]
for m in methods:
    optim_res = optimize.minimize(fun, x0, method=m)
    print("---\nMethod:{}\n".format(m),optim_res, "\n")
Explanation: On voit que le minimum se trouve près de $[-\frac{1}{2}, 0]$. On va utiliser ce point pour initialiser l'optimisation.
On va tester différentes méthodes et comparer les sorties obtenues.
End of explanation
for i in range(4):
    optim_res = optimize.minimize(fun, x0, method='BFGS')
    print("---\nMethod:{} - Test:{}\n".format(m,i),optim_res, "\n")
Explanation: On trouve un minimum à $-0.4052$ en $[-0.669, 0.000]$ pour toutes les méthodes qui convergent. Notez le message de sortie de 'CG' qui signifie que le gradient ne varie plus assez. Personnellement, je ne trouve pas ce message de sortie très clair. Le point trouvé est bien l'optimum cherché pourtant. Notez aussi le nombre d'évaluations de la fonction (nfev) pour chaque méthode, et le nombre d'évaluation de gradient (njev) pour les méthodes qui reposent sur un calcul de gradient.
Remarquez aussi que si on relance Anneal plusieurs fois, on n'est pas assuré d'obtenir la même solution, puisqu'il s'agit d'une métaheuristique.
End of explanation
for m in methods:
    print("Method:{}:".format(m))
    %timeit optim_res = optimize.minimize(fun, x0, method=m)
    print('############')
Explanation: On va évaluer le temps de calcul nécessaire à chaque méthode.
End of explanation
def shifted_scaled_bowlpeak(x,a,b,c):
    return (x[0]-a)*np.exp(-((x[0]-a)**2+(x[1]-b)**2))+((x[0]-a)**2+(x[0]-b)**2)/c
a = 2
b = 3
c = 10
optim_res = optimize.minimize(shifted_scaled_bowlpeak, x0, args=(a,b,c), method='BFGS')
print(optim_res)
print('#######')
optim_res = optimize.minimize(lambda x:shifted_scaled_bowlpeak(x,a,b,c), x0, method='BFGS')
print(optim_res)
Explanation: On peut aussi fournir des arguments supplémentaires à la fonction qu'on optimise. Par exemple, les données lorsque vous maximisez une log-vraissemblance. En voici un exemple: on considère une version rescaled de la fonction bowl_peak. Vous pourriez aussi utiliser une lambda fonction.
End of explanation
<END_TASK>