Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

A new filter #382

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Conversation

richardalligier
Copy link

Adding FilterConsistency, a filter based on finding the longest sequence of consistent points.




def distance(lat1,lon1,lat2,lon2):
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Il y a une dépendance https://github.com/open-aviation/pitot
Dedans tu as from pitot.geodesy import distance
Pour les unités, tu peux utiliser le décorateur @impunity avec les unités en annotations.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

J'avais vu et j'avais essayé le distance de pitot, mais j'ai une exception qui est lié au fait que lat1 et lat2 n'ont pas même shape. Sur un exemple lat1.shape=(200,848) tandis que lat2.shape=(848,).

On peut s'en sortir en broadcastant avec np.broadcast_to avant l'appel à pitot.geodesy.distance mais c'est long en temps de calcul: 0.135s. Par contre en utilisant ma distance: 0.043s.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Après ta distance, c'est la formule de Haversine en sphérique.
Ton dx/dy, c'est une formule en ellipsoïdal.

Ça expliquerait le delta de temps, et ça reboucle sur ma question de faire des calculs en euclidien après projection en local (là au moins les calculs sont plus simples/rapides).

Si la précision est importante, pitot fait appel à pyproj derrière (c'est juste des trucs malins de typage + unités en plus), et c'est quand même assez optimal (et maintenu par des gens plus intelligents que moi ^^)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ouai, ça pourrait être une idée de projeter au tout début mais j'ai un peu peur pour les longues trajectoires. Et puis je sais pas c'est quoi le pire entre l'erreur de distance sur des projections et l'erreur d'Haversine.

Sur ma trajectoire test, l'écart max est de 0.17% en relatif entre Haversine et geodesy.distance. J'ai pas trouvé d'équivalent WGS84 sympa à Haversine.

Pour le dx,dy, j'ai la formule ellipsoid donc je l'utilise.

Bref, à réfléchir...

res = np.empty((horizon,v.shape[0]))
res.fill(np.nan)
for i in range(horizon):
res[i,:v.shape[0]-i] = v[i:]
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

C'est pas long cette boucle? np.pad fait pas le boulot plus rapidement?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

En fait c'est pas un padding que je fais, c'est surtout un shift qui me permet d'avoir res[h,i] = v[i+h].
C'est intéressant, car ensuite si je fais diff=lag(horizon,v)-v, avec le broadcasting je vais avoir diff[h,i] = v[i+h]-v[i] pour tout h<horizon.

Pour éviter les boucles, il y a moyen de faire de la "magie noire" avec les np.lib.stride_tricks. Une façon de faire effectivement plus rapide est: np.lib.stride_tricks.sliding_window_view. Une façon encore plus rapide est de faire avec np.lib.stride_tricks.as_strided mais celle-ci est accompagné d'un "This function has to be used with extreme care, see notes.". Ah, et au final j'utilise np.pad pour mettre des np.nan en fin d'array avant de faire des stride_tricks.

Sur ma trajectoire exemple:
_l'ancienne version prend 0.28ms par appel de lag;
_ la version avec sliding_window_view prend 0.15ms par appel;
_ la version avec as_strided prend 0.10ms par appel.

Au final ça change pas la guerre comparé au 40ms total de l'algo. Mais bon, c'est toujours bon à prendre, du coup je pars sur np.lib.stride_tricks.sliding_window_view !

d = a - b
return d + 2 * np.pi * ((d<-np.pi).astype(float)-(d>=np.pi).astype(float))

def dxdy_from_dlat_dlon(lat_rad, lon_rad, dlat, dlon):
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Est-ce que ça vaut pas le coup de passer ça dans pitot/geodesy?
C'est une formule qui ne marche que localement? Est-ce que c'est plus efficace que de faire ce que je fais souvent d'habitude, i.e. calculer x,y sur une projection locale et bosser en euclidien?

@xoolive
Copy link
Owner

xoolive commented Nov 29, 2023

Merci!

Je mets quelques commentaires (résumé du ligne à ligne au dessus) pour ne pas oublier ici :

  • est-ce que tu pourrais lancer poetry run ruff format src avant de recommitter;
  • je ferai les annotations de type plus tard si tu ne sais pas trop faire ;
  • les dépendances https://github.com/achevrot/impunity et https://github.com/open-aviation/pitot sont déjà là, donc autant réutiliser ;
  • est-ce que tu aurais des exemples pour illustrer ? (qui pourraient faire à la fois tests et documentations)

Et pour graph-tool: https://git.skewed.de/count0/graph-tool/-/wikis/installation-instructions#native-installation

@xoolive
Copy link
Owner

xoolive commented Nov 29, 2023

J'ai aussi corrigé le test qui marchait plus sur le repo d'origine; tu devrais faire un rebase.

@xoolive
Copy link
Owner

xoolive commented Apr 13, 2024

autre bibliothèque à tester aussi: https://www.rustworkx.org/apiref/rustworkx.dag_longest_path.html

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants