Comment optimiser ce code ?

Bonjour à tous,

je suis un novice dans la programmation python et j'aurai aimé savoir comment je pourrais faire pour optimiser le code suivant :

Code:

---

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

import psyco psyco.full()   obs = [0, 0, 1, 1, 0, 1] pred = [0.9874562164, 0.23565464, 0.756586787, 0.894542311, 0.1212121, 0.9953121234]   def sum_sens_spec(obs, pred): vp = 0 vn = 0 fp = 0 fn = 0   b = xrange(len(obs))   for j in b: if obs[j]: if pred[j]: vp +=1 else: fn +=1 else: if pred[j]: fp +=1 else: vn +=1   sens = vp/(vp+fn) spec = vn/(vn+fp) res1 = sens+spec return res1   cuts = [u/n for u in xrange(n+1)] res = []   for u in cuts: pred2 = [i>=u for i in pred] res.append(sum_sens_spec(obs, pred2))   m1 = max(res)   u = 0 while res[u] < m1: u +=1

---

L'objectif de ce code est de trouver la valeur comprise entre 0 et 1 qui permet de maximiser la somme de spec et sens. Ici obs et pred prennent que 6 valeurs mais dans la réalité ils en prennent plus de 1200 en moyenne.

Toutes idées pour optimiser le temps de calcul serait bienvenue.

Merci

Edit : je bosse avec python 2.4.2

Il fonctionne ce code ??

Code:

---

cuts = [u/n for u in xrange(n+1)]

---

En dehors du fait que n n'est pas défini, en Python 2.x, / sera ici une division entière et donnera donc 0 pour tous les éléments de cuts sauf le dernier, et je doute que ce soit l'intention.

Code:

---

1 2 3 4 5 6

>>> n = 6 >>> [u/n for u in xrange(n+1)] [0, 0, 0, 0, 0, 0, 1] >>> [float(u)/n for u in xrange(n+1)] [0.0, 0.16666666666666666, 0.3333333333333333, 0.5, 0.6666666666666666, 0.8333333333333334, 1.0] # ce serait pas plutôt ça l'intention ?

---

Même remarque pour les autres divisions, c'est pas plutôt une division flottante que tu veux ?

Ou alors tu as fait un from __future__ import division précédemment ? (Je ne sais même pas si cela existe en Python 2.4).

Ce sera difficile de t'aider si on n'a pas un code correct et complet...

Bon, en supposant que le code est correct modulo mes remarques ci-dessus, voilà ce à quoi je suis arrivé:

Il s'agit principalement de tirer parti du fait que "obs" ne change pas à chaque appel de sum_sens_spec.

1) division par vp+fn: vp+fn est égal au nombre de 1 dans obs, qui est constant, on peut donc le précalculer (j'ai appelé la variable fnvp pour rappeler d'où ça vient). Idem pour vn+fp, qui est le nombre de 0 dans obs.

2) le test if obs[j] peut être évité si on partitionne à l'avance pred en 2 morceaux (pred0 et pred1), qui correspondent aux valeurs de pred pour lesquelles les valeurs correspondantes de obs sont respectivement 0 et 1.

Dès lors, obs n'est plus utile dans sum_sens_spec.

3) Les boucles du genre:

Code:

---

1 2 3 4 5

vp = 0 b = xrange(len(obs)) for j in b: if pred1[j]: vp +=1

---

peuvent être remplacées par:

Code:

---

vp = sum(pred1)

---

vu que pred1 est une liste de booléen (après avoir été filtrée par le cut). Voire même, en y intégrant directement le cut:

Code:

---

vp = sum([i>=u for i in pred1]) # avec u = valeur courante du cut

---

Au final, sum_sens_spec devient tellement simple que l'on peut se passer d'en faire une fonction séparée et on gagne ainsi un appel de fonction.

Au final cela donne ceci (j'ai généré aléatoirement des inputs pour tester):

Code:

---

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 29 30 31 32 33 34 35

from __future__ import division n = 2000 # nombre de cuts m = 2000 # taille de obs et de pred   #obs = [0, 0, 1, 1, 0, 1] #pred = [0.9874562164, 0.23565464, 0.756586787, 0.894542311, 0.1212121, 0.9953121234]   # pour simuler des inputs import random random.seed(0) obs = [random.randint(0,1) for _ in xrange(m)] pred = [random.random() for _ in xrange(m)]     # précalculer le nombre de 0 et de 1 dans obs fnvp = sum(obs) fpvn = len(obs) - fnvp   # partitionner pred en 2 morceaux: # pred0 qui contient les valeur de pred pour lesquelles obs=0 # pred1 qui contient les valeur de pred pour lesquelles obs=1 pred = zip(obs, pred) # decorate pred.sort() # sort pred = [p for o,p in pred] # undecorate pred0, pred1 = pred[:fpvn], pred[fpvn:] # partitionnement   cuts = [u/n for u in xrange(n+1)]   # calcul de res (nettement simplifié par les manipulations précédentes) res = [sum([i>=u for i in pred1]) / fnvp + sum([i<u for i in pred0]) / fpvn for u in cuts]   m1 = max(res) u = res.index(m1) # équivalent de la boucle while   print m1, u

---

Il y a peut-être moyen d'optimiser un peu plus en Python pur, sinon numpy devrait permettre d'aller plus loin...

Bonjour,

Effectivement j'ai oublié quelques lignes de codes. J'avais dans mon script spécifié from __future__ import division ainsi que n qui vaut en réalité 20000.

En tout cas merci pour tous ces conseils je vais me plongé dans la lecture de ton code de ce pas.

Encore merci, j'ai appris plein de choses.

Malheureusement après un test à "données réelles" le code que tu me proposes tournes plus lentement que le code que j'avais utilisé.

Merci encore.

Euh tu as rajouté

Code:

---

1 2	import psyco psyco.full()

---

à son code avant de faire la comparaison ?
Sinon ça n'a pas de sens.

Oui j'avais bien rajouté ces deux lignes. Mes comparaisons ont été menées dans les mêmes conditions.

Pourtant cela tournait 2 à 3 fois plus vite dans mes tests (en Python 2.7, donc sans psyco, vu qu'il n'a pas encore été porté...).

Mais je ne suis pas tellement étonné; optimiser pour CPython et optimiser pour psyco sont deux choses différentes. psyco est surtout efficace pour optimiser les appels de fonctions et les calculs arithmétiques, et c'est principalement ce que j'ai retiré. Optimiser pour CPython (sans psyco), c'est souvent essayer d'exploiter au maximum les builtins et fonctions de bibliothèques qui sont déjà codées en C, psyco n'a aucun effet là-dessus...

tototode, pourrais-tu décrire brièvement à quoi correspond cet algorithme ? Que représente les flottants dans pred ? des probabilités ? Et obs, je suppose qu'il s'agit d'observations ? Qu'est-ce qui est calculé avec cet algorithme ?

Les optimisations que j'ai faites, c'est du "bas-niveau", des transformations de programmes à "sémantique constante", mais savoir ce qu'on calcule ça peut aider pour trouver une meilleur approche...

Une percée algorithmique:

Le nombre de valeurs de pred0 qui sont < à la valeur courante de cut est une fonction monotone de cut. Donc, quand la valeur de cut croit, il est inutile de recalculer les éléments qui ont déjà été calculés. Si, de plus, les valeurs de pred0 sont triées (ce qui est déjà le cas dans mon code ci-dessus), on peut se contenter d'un seul parcours de chaque tableau (pred0 et pred1). Pour pred1 le parcours doit être fait dans l'autre sens (décroissant).

Cela donne ceci:

Code:

---

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53

from __future__ import division n = 20000 # nombre de cuts m = 1200 # taille de obs et de pred   #obs = [0, 0, 1, 1, 0, 1] #pred = [0.9874562164, 0.23565464, 0.756586787, 0.894542311, 0.1212121, 0.9953121234]   # pour simuler des inputs import random random.seed(0) obs = [random.randint(0,1) for _ in xrange(m)] pred = [random.random() for _ in xrange(m)]   # précalculer le nombre de 0 et de 1 dans obs fnvp = sum(obs) fpvn = len(obs) - fnvp   # partitionner pred en 2 morceaux: # pred0 qui contient les valeur de pred pour lesquelles obs=0 # pred1 qui contient les valeur de pred pour lesquelles obs=1 pred = zip(obs, pred) # decorate pred.sort() # sort pred = [p for _,p in pred] # undecorate pred0, pred1 = pred[:fpvn], pred[fpvn:] # partitionnement   cuts = [u/n for u in xrange(n+1)]   # calcul de res pred0s = [] i = 0 l = len(pred0) for u in cuts: while i < l and pred0[i] < u: i += 1 pred0s.append(i)   pred1.reverse() # on veut pred1 décroissant cuts.reverse() # et cuts décroissant pred1s = [] i = 0 l = len(pred1) for u in cuts: while i < l and pred1[i] >= u: i += 1 pred1s.append(i) pred1s.reverse() # mettre pred1s en ordre décroissant   res = [p1 / fnvp + p0 / fpvn for p0, p1 in zip(pred0s, pred1s)]   m1 = max(res) u = res.index(m1)   print m1, u

---

Dans mon test, le temps de calcul est divisé par 100 environ (avec 20000 cuts et 1200 inputs). On passe d'un algo qui est O(m*n) à O(m log m + n) si je ne m'abuse.

Dans le calcul de pred1s, les appels à reverse() pourraient être éliminés en parcourant cuts avec une boucle while (au lieu d'un for) et des indices décroissants. Le gain serait certainement nul en Python pur mais peut-être avec psyco il serait mesurable...

re,

en fait les obs sont des valeurs de présence absence d'une espèce dans un site donné et les preds sont les prédictions issues d'une régression logistique.

L'objectif de ce code est de trouver la probabilité "seuil" qui maximise la somme de la sensibilité (taux de présence correctement prédite) et de la spécificité (taux d'absence bien prédite).

En tout cas merci de l'intérêt que tu portes à cet algorithme. Je vais tester ceci de suite.

Le dernier code est beaucoup plus efficace et le gain de temps est relativement important. Encore merci à toi, à la fois pour avoir pris le temps de répondre, pour le code qui fonctionne, mais aussi pour tout ce que tu m'as appris.

Merci pour l'explication. Au moins je sais sur quoi j'ai travaillé ;)

[EDIT]
Avec le recul, je me demande à quoi cela sert de tester 20000 seuils différents alors qu'il n'y a qu'environ 1200 inputs. Il y aura beaucoup de ces valeurs testées qui donneront la même sensibilité et la même spécificité. Par exemple, après avoir exécuté le code précédent:

Code:

---

1 2 3 4

>>> len(res) 20001 >>> len(set(res)) # nombre de valeurs distinctes 1129

---

Pourquoi ne pas utiliser les valeurs de pred comme seuils possibles ? On arrivera au même résultat (le maximum sensibilité+spécifité trouvé sera le même*; le seuil quant à lui sera toujours une valeur présente dans pred).
Il suffit de remplacer:

Code:

---

cuts = [i/n for i in xrange(n+1)]

---

par

Code:

---

1 2	cuts = pred cuts.sort()

---

Non ?

* En théorie, il pourrait même être meilleur

Au départ j'étais parti sur 20000 seuils parce que j'ai vu graphiquement que la somme de la sensibilité et de la spécificité présentait de nombreuses micro-variations et qu'en deçà de ce seuil on pouvait avoir des estimations légèrement erronées. Au départ j'avais utilisé un algorithme de minimisation mais le seuil estimé bien que proche des seuils que je calcule avec ton code, ne donnait pas la somme maximale de la sensibilité et de la spécificité.
Il est important de disposer d'une estimation fiable de ce seuil car ces résultats sont intégrés dans une publication dans laquelle on dit que le seuil est calculer de sorte à maximiser la somme de la spécificité et de la sensibilité. Donc autant l'estimer avec la plus grande fiabilité possible.

J'avais aussi pensé à d'autres pistes pour optimiser le temps de calcul comme par exemple, commencer par faire ce calcul pour un écart de 1/100, trouver le seuil qui maximise la somme et ensuite regarder plus finement ce qui se passe autour de ce seuils (dans un intervalle de + ou - 0.01 par exemple).

Maintenant le temps de calcul est super satisfaisant donc je remiserais ces améliorations pour plus tard.

Merci

Se limiter aux valeurs de prédictions calculées ne diminue en rien la fiabilité du résultat, au contraire. Comme le nombre de point de données est fini, la fonction résultat est une fonction en escalier sur l'intervalle [0,1]. En tout cas, cette approche est cohérente avec le code original; on peut imaginer interpoler les valeurs entre ces points de données, mais ce serait pousser le bouchon un peu loin. Il y a donc un intervalle de valeurs pour lesquelles la fonction atteint son maximum. Cet intervalle est forcément délimité par deux points de données, c'est à dire deux valeurs présentes dans les prédictions. Utiliser les points de données comme seuils possibles permet de déterminer cet intervalle avec précision; cette solution est optimale.

En utilisant des seuils régulièrement espacés, il y a toujours la possibilité de "sauter au-dessus" de cet intervalle, s'il est plus court que l'espacement entre les seuils.

Dans mon post précédent, j'avais indiqué comme tirer parti de cela avec une modification minimale du code, mais en fait cela permet de simplifier nettement le code, et de le rendre encore plus rapide.

Cette (dernière, je pense) version affiche l'intervalle de valeurs parmi lequel choisir le seuil pour maximiser sensibilité+spécificité.

Code:

---

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

from __future__ import division m = 1200   import random random.seed(0) pred = [random.random() for _ in xrange(m)] obs = [random.random() <= p for p in pred]   nb_present = sum(obs) nb_absent = len(obs) - nb_present   pred_obs = zip(pred, obs) pred_obs.sort() pred, obs = zip(*pred_obs)   res = [] val = 0 for o in obs: res.append(val) val += [nb_present, -nb_absent][o] m1 = max(res) u = res.index(m1) m1 = m1 / (nb_absent * nb_present) + 1.0   print "max=%f cut=]%f,%f]" % (m1, pred[u-1] if u>0 else 0, pred[u])

---

J'ai fait un test en utilisant 1000 jeux de données générés aléatoirement comme ci-dessus (j'ai changé la génération aléatoire par rapport à mes premiers posts pour qu'il y ait un minimum de corrélation entre les observations et les prédictions), en mettant en concurrence les deux façons de faire (seuils régulièrement espacés (algo1) vs intervalle de seuils déterminés par les prédictions (algo2)).

Dans 965 cas, le seuil déterminé par l'algo1 se trouve dans l'intervalle déterminé par l'algo2 (et le maximum calculé est le même).

Dans 31 cas, le maximum de l'algo2 est meilleur que celui de l'algo1 (et le seuil n'est donc pas dans l'intervalle). Cela montre la supériorité de l'algo2.

Dans les 4 derniers cas, les maxima des 2 algos sont identiques mais le seuil de l'algo1 n'est pas dans l'intervalle de l'algo2. Ces cas sont imputables à la précision limitée des floats; même dans ce cas je fais plus confiance à l'algo 2 car tous les calculs sont effectués sur des entiers, le max n'est converti en flottant qu'à la fin pour comparaison.


Bon voilà. Cela dit, tu n'es pas obligé de me croire ni de changer ton code; mais je ne voulais pas te laisser sur l'impression que cette approche est moins fiable que la précédente ;)

Je partage ton avis concernant l'utilisation des probabilités prédites plutôt que des seuils fixés à un pas de temps donné (1/n ici). En plus de simplifier le calcul, cette solution semble plus élégante et surtout plus en lien avec l'objectif initial et en adéquation avec les données.
Je vais comparé sur mes données les résultats issus des deux algos pour voir les différences.

Encore merci pour tout.

Citation:

---

Code:

---

print "max=%f cut=]%f,%f]" % (m1, pred[u-1] if u>0 else 0, pred[u])

---

---

Je me permets de revenir sur ton dernier code. Je ne comprends pas très bien ce que tu veux faire passer avec cette ligne de code et quand utiliser pred[u-1] plutôt que pred[u].

Encore merci.

Citation:

---

Envoyé par tototode

Je me permets de revenir sur ton dernier code. Je ne comprends pas très bien ce que tu veux faire passer avec cette ligne de code et quand utiliser pred[u-1] plutôt que pred[u].

---

La réponse est un intervalle pred[u-1], pred[u]. N'importe quel seuil choisi dans cet intervalle donnera la valeur optimale pour la mesure (sens+spec). Ne fais pas trop attention au fait que cet intervalle est ouvert d'un côté et fermé de l'autre, cela n'est le fait que d'un "tuning" pour qu'il respecte la spécification donnée par ton code original (le choix de >= au lieu de > dans [i>=u for i in pred]). Il vaut sans doute mieux éviter les extrémités et choisir un seuil entre ces deux valeurs, par exemple (pred[u-1]+pred[u])/2.