# TD3 : Numpy
## Partie 1 : Approche naïve du calcul scientifique en Python
Le fichier `basic_func.py` vous est fourni et il contient une fonction ainsi qu'un `main` :
```python
import numpy as np
import time
def average(a):
"""
Compute the average value of ``a``
Parameters
----------
a : ndarray
Input array.
Returns
-------
avg : int
The average value of ``a``
"""
average = 0
for i in range(a.shape[0]):
average += a[i]
return average/a.shape[0]
if __name__ == "__main__":
rand_array_1 = np.random.randint(0, 10000, 1000000)
t0 = time.time()
average_value = average(rand_array_1)
print(time.time() - t0)
t0 = time.time()
average_value = np.mean(rand_array_1)
print(time.time() - t0)
```
Questions :
- Que fait la fonction `average` ?
- Que font les lignes du `main` ?
- Quels résultats obtient-on ?
- Qu'en déduire ?
## Partie 2 : Généralités sur les tableaux `numpy`
### Exercice 1
Dans la suite de cet exercice, le tableau de données `foo` sera manipulé. Il est créé comme suit :
```python
foo = np.random.randint(
0,
10,
(np.random.randint(50, 60), np.random.randint(42, 69))
)
```
1. Que fait cette ligne ?
**Réponse :** Cette ligne crée une matrice d'entiers entre 0 et 10 de taille aléatoire : entre 50 et 59 lignes et entre 42 et 68 colonnes.
2. Comment connaître le **nombre de lignes** d'un tableau ?
**Réponse :** L'attribut `foo.shape` est un tuple (nb de lignes, nb de colonnes). `foo.shape[0]` renvoie donc le nombre de lignes.
3. Comment connaître le **nombre de colonnes** d'un tableau ?
**Réponse :** `foo.shape[1]` renvoie donc le nombre de colonnes.
4. Comment connaître le **nombre d'éléments** d'un tableau ?
**Réponse :**`foo.size` renvoie la taille d'un tableau (son nombre d'éléments)
5. Comment connaître le **type des éléments** d'un tableau ? Et le changer ?
**Réponse :** L'attribut `foo.dtype` renvoie le type des éléments de `foo`. `foo.astype('float64')` ou `foo.astype(np.float)` permet de convertir les éléments vers un autre type.
6. Que renvoie `foo > 8` ?
**Réponse :** Un tableau de booléens qui est obtenu en appliquant la condition à chaque élément.
7. Comment **extraire les coordonnées** des éléments strictement inférieurs à 5 ?
**Réponse :** La fonction `where(condition)` permet d'extraire les coordonnées des éléments qui vérifient cette condition.
8. Comment **multiplier ces éléments par 2** ?
**Réponse :** On peut accéder aux éléments étant donné une liste de coordonnées `foo[np.where(foo < 5)]`. On peut donc les multiplier en faisant `foo[np.where(foo > 5)] *= 2`.
9. Que font les lignes suivantes ?
```python
bar = foo[:, 0:13:2]
bar_copy = foo[:, 0:13:2].copy()
bar_float = bar.astype(np.float)
baz = foo[4:14,(1,6,11,26)]
```
**Réponses :**
- *ligne 1 :* Récupération de toutes les lignes et une colonne sur deux de 0 à 13.
- *ligne 2 :* Récupération des mêmes informations mais force la copie des données (modifier `bar` ne modifiera pas `foo`).
- *ligne 3 :* Conversion en flottants.
- *ligne 4 :* Récupération de la ligne 4 à 14 et les colonnes 1, 6, 11 et 26.
### Exercice 2
1. Calculer l'hypothénuse de chaque triangle en utilisant le théorème de Pythagore.
**Solution :**
```python
hypothenuse = np.sqrt(short_sides[:,1]**2+short_sides[:,0]**2)
```
2. Créer `triangles` qui rassemble `short_sides` et les hypothénuses.
**Solution :**
```python
triangles = np.column_stack((short_sides, hypothenuse))
```
3. Calculer un des deux angles non droits via **l'arc cosinus**, **l'arc sinus** et **l'arc tangente**.
**Solution :**
```python
acos_tr = np.arccos(triangles[:,0]/triangles[:,2])
asin_tr = np.arcsin(triangles[:,1]/triangles[:,2])
atan_tr = np.arctan2(triangles[:,1], triangles[:,0])
```
4. Vérifier si les valeurs trouvées sont **toutes** égales ?
**Solution :**
```python
np.all(acos_tr==asin_tr)
np.all(acos_tr==atan_tr)
```
5. Vérifier si **quelques** valeurs trouvées sont égales ?
**Solution :**
```python
np.any(acos_tr==asin_tr)
np.any(acos_tr==atan_tr)
```
6. Vérifier si les valeurs trouvées sont **à peu près** égales ?
**Solution :**
```python
np.allclose(acos_tr, asin_tr)
np.allclose(acos_tr, atan_tr)
```
7. Que signifie **à peu près** ?
**Réponse :** Par défaut, la tolérance relative est de 1e-5 et la tolérance absolue de 1e-8 mais peut être modifiée en utilisant la fonction comme suit :
```python
numpy.allclose(acos_tr, asin_tr, rtol=1e-8, atol=1e-10)
```
La fonction vérifie la condition suivante :
```python
numpy.allclose(a, b, rtol=1e-5, atol=1e-8)
# absolute(a - b) <= (atol + rtol * absolute(b))
```
## Partie 3 : Millions de mesures de coquilles
### Préliminaire : Chargement de données
Pour charger vos données :
- Utilisez tout d'abord `np.loadtxt` pour charger un fichier texte dans un tableau numpy.
- Puis `np.save` permet de sauver les données déjà chargées dans tableau dans un format intermédiaire.
- Enfin, vous pouvez charger vos données avec `np.load` pour gagner du temps par rapport à la première étape.
```python
# /!\ A NE FAIRE QU'UNE FOIS
data=np.loadtxt("coquilles.txt")
# /!\ A NE FAIRE QU'UNE FOIS
np.save("coquilles.npy", data)
# A répéter au lancement du programme
data = np.load("coquilles.npy")
```
Le tableau contient les données formattées comme suit :
- Colonne 1 : Identifiant de la coquille (0 à 13400)
- Colonne 2 : Époque d'enregistrement (0 à 196)
- Colonne 3 : Latitude de la coquille
- Colonne 4 : Longitude de la coquille
- Colonne 5 : Stade d'évolution (1 ou 5)
Le début du fichier ressemble à :
```python
[ 0. 0. -0.19542887 49.89834643 1. ]
[ 0. 1. -0.17596363 49.91657752 1. ]
...
```
### Migration de coquilles
---
Pour tracer la trajectoire de la 42e coquille :
- Créer un masque d'indices qui permette de déterminer les positions de la 42e coquille
- Extraire la lattitude et la longitude
- Afficher avec `matplotlib`
**Solution :**
```python
# Extraction des données
mask_42 = data[:, 0] == 41
longi_42 = data[mask_42][:, 2]
lat_42 = data[mask_42][:, 3]
# Affichage
plt.figure()
plt.plot(longi_42, lat_42)
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.title("42th scallop trajectory")
```
---
Pour retrouver les coquilles en stage 5 avant le 60e enregistrement il faut recouper deux conditions. On peut recouper deux valeurs booléennes avec un `*` (opération `AND`) ou `+` (opération `OR`).
```python
a = np.array([True, True, True])
b = np.array([True, False, False])
print(a*b)
>>> np.array([True, False, False])
print(a+b)
>>> np.array([True, True, True])
```
**Solution :**
```python
mask_early = (data[:, -1] == 5) * (data[:, 1] < 60)
early = np.unique(data[mask_early][:, 0])
print(early)
```
---
Le point de départ de chaque coquille correspond aux coordonnées à l'étape 0.
**Solution :**
```python
start = data[data[:, 1] == 0][:, 2:4]
```
---
Pour trouver le changement de stage 1 au stage 5 on peut :
- Extraire les stages des coquilles
- Copier en décalant d'un indice les stages (`np.roll`)
- Chercher quand la différence des deux vaut `-4` (donc un changement de stage)
```python
x = np.arange(10)
np.roll(x, 2)
>>> array([8, 9, 0, 1, 2, 3, 4, 5, 6, 7])
np.roll(x, -2)
>>> array([2, 3, 4, 5, 6, 7, 8, 9, 0, 1])
```
Pour pouvoir les trier, la fonction `argsort()` renvoie les indices de manière à ce que le tableau résultant soit trié.
```python
x = np.array([3, 1, 2])
np.argsort(x)
>>> array([1, 2, 0])
```
**Solution :**
```python
stage_col = data[:, -1]
stage_col_roll = np.roll(stage_col, 1)
mask_stage = (stage_col_roll - stage_col) == -4 # !!! parenthèses !!!
stage_change = data[mask_stage]
# trier les coquilles en fonction de leur date de changement d'étape
stage_sorted = stage_change[stage_change[:, 1].argsort()]
print(stage_sorted[:, 0].astype(np.int))
```
---
###### tags: `python`
Sign in with Wallet
Connect another wallet