Algorithme pour trouver la masse agrégée de structures semblables à des «barres granola»?

Je suis un chercheur en sciences planétaires et un projet sur lequel je travaille est la simulation à N corps des anneaux de Saturne. Le but de cette étude particulière est de regarder les particules s'agglutiner sous leur propre gravité et de mesurer la masse agrégée des mottes en fonction de la vitesse moyenne de toutes les particules de la cellule. Nous essayons de savoir si cela peut expliquer certaines observations faites par le vaisseau spatial Cassini pendant le solstice d'été saturnien lorsque de grandes structures ont été vues projetant des ombres sur les anneaux presque bordés. Voici une capture d'écran de ce à quoi ressemble un pas de temps donné. (Chaque particule a un diamètre de 2 m et la cellule de simulation elle-même mesure environ 700 m de diamètre.)

Le code que j'utilise déjà crache la vitesse moyenne à chaque pas de temps. Ce que je dois faire, c'est trouver un moyen de déterminer la masse de particules dans les mottes et NON les particules parasites entre elles. Je connais la position, la masse, la taille, etc. de chaque particule, mais je ne sais pas facilement que, disons, les particules 30 000-40 000 ainsi que 102 000-105 000 constituent un brin qui, à l'œil humain, est évident.

Ainsi, l'algorithme que je dois écrire devrait être un code avec aussi peu de paramètres saisis par l'utilisateur que possible (pour la réplicabilité et l'objectivité) qui passerait par toutes les positions des particules, déterminer quelles particules appartiennent aux amas, puis calculer la Masse. Ce serait formidable s'il pouvait le faire pour "chaque" touffe / brin par opposition à tout sur la cellule, mais je ne pense pas en avoir vraiment besoin pour les séparer.

La seule chose à laquelle je pensais était de faire une sorte de calcul de la distance N ² où je calculerais la distance entre chaque particule et si, disons, les 100 particules les plus proches se trouvaient dans une certaine distance, alors cette particule serait considérée comme faisant partie d'un grappe. Mais cela semble assez bâclé et j'espérais que vous, les gens et les programmeurs CS, connaissiez une solution plus élégante?

Modifié avec ma solution: ce que j'ai fait, c'est d'adopter une sorte d'approche de voisin / cluster le plus proche et de commencer par l'implémentation de N ² rapide et sale . Donc, prenez chaque particule, calculez la distance par rapport à toutes les autres particules, et le seuil pour dans un cluster ou non était de savoir s'il y avait N particules dans la distance d (deux paramètres qui doivent être définis a priori , malheureusement, mais comme certains l'ont dit réponses / commentaires, je n'allais pas m'en tirer sans en avoir).

Je l'ai ensuite accéléré en ne triant pas les distances, mais en faisant simplement une recherche d' ordre N et en incrémentant un compteur pour les particules dans d , et cela a accéléré les choses d'un facteur 6. Puis j'ai ajouté un "arbre de programmeur stupide" (parce que je sais presque rien sur les codes d'arbre). Je divise la cellule de simulation en un nombre défini de grilles (meilleurs résultats lorsque la taille de la grille ≈7 d ) où la grille principale s'aligne avec la cellule, une grille est décalée de moitié en x et y , et les deux autres sont décalées de 1/4 en ± x et ± y . Le code divise ensuite les particules en grilles, puis chaque particule N n'a qu'à avoir des distances calculées par rapport aux autres particules de cette cellule.

Théoriquement, s'il s'agissait d'un véritable arbre, je devrais obtenir la commande N * log ( N ) par opposition à N ² vitesses. Je suis arrivé quelque part entre les deux, où pour un sous-ensemble de 50 000 particules, j'ai obtenu une augmentation de 17x de la vitesse, et pour une cellule de 150 000 particules, j'ai eu une augmentation de 38x de la vitesse. 12 secondes pour la première, 53 secondes pour la seconde, 460 secondes pour une cellule de 500 000 particules. Ce sont des vitesses comparables au temps que prend le code pour exécuter le pas de temps de simulation 1, c'est donc raisonnable à ce stade. Oh - et il est entièrement fileté, il faudra donc autant de processeurs que je peux y jeter.

Ma première suggestion est de diviser votre problème en deux problèmes: d'abord, déterminez ce que vous voulez, puis trouvez comment obtenir efficacement ce que vous voulez. Vous ne pouvez pas obtenir efficacement quelque chose que vous n'avez pas encore défini. Je vais mettre quelques idées dans cette réponse qui peuvent vous aider à trouver cette définition. Je vous suggère de faire une mise en œuvre inefficace des idées que vous aimez d'abord, de l'appliquer à quelques ensembles de données pas trop grands, d'évaluer les résultats à la main, d'adapter votre définition et de répéter (éventuellement en posant une autre question ici), jusqu'à ce que vous soyez satisfait de votre définition. Après cela, je vous suggère de poser une autre question sur la façon de calculer efficacement le résultat de votre définition (si vous avez encore besoin d'aide).

Voyons donc ce qui correspondrait à notre idée intuitive d'un «brin». Vos brins semblent être constitués de points à peu près uniformément distribués, mais vous devriez le vérifier en faisant une image agrandie (du jeu de données d'origine) - la résolution de votre image est trop faible pour dire avec certitude que les points sont vraiment à peu près uniformément distribués . Je suppose qu'ils sont pour cette réponse.

Une première idée peut être de regarder le plus proche voisin de chaque point. Choisissons un point X, appelons son voisin le plus proche Y et définissons D comme la distance entre X et Y. Nous regardons ensuite le cercle C autour de X de rayon D * A, où A est un paramètre de réglage, disons A = 3. Si X fait partie d'un brin, nous nous attendons à ce que pour chaque point Z dans C, la distance de Z à son plus proche voisin W soit à peu près la même que D. Si elle est significativement plus courte, disons plus que A (ou peut-être un autre paramètre B) alors X est apparemment proche de points beaucoup plus proches les uns des autres que de X, donc X ne fait probablement pas partie d'un brin.

Ce critère n'est cependant pas complet. Il ne donne qu'un critère pour détecter une «frontière» entre des zones denses à points et des zones moins denses à points. Nous devons encore regrouper les points en brins.

Il y a une fonctionnalité dans votre image qui montre que ce n'est pas simple. Dans le coin inférieur droit de votre image, il y a une zone relativement grande avec beaucoup de points parasites. Ces points parasites sont eux-mêmes à peu près uniformément distribués, donc si nous supprimions tous les points du brin qui l'entoure (et tous les autres points), nous nous attendrions à ce que tout algorithme de détection de brin marque cet ensemble de points parasites comme un brin! Nous devons donc être prudents lors de la création de nos clusters.

Une idée peut être de faire ce qui suit. Nous allons faire un graphique sur ces points, où les sommets sont les points et les bords signifient que deux points ont une densité similaire. Pour chaque point, nous vérifions le critère ci-dessus. S'il vérifie, nous connectons X avec un bord à tous les points de C. S'il ne récupère pas, nous n'ajoutons aucun bord et marquons X comme «parasite». Après avoir fait cela pour chaque point, nous considérons l'ensemble des composants connectés. Ceux-ci doivent consister en un seul (dans le cas de votre image, mais d'autres ensembles de données peuvent avoir plusieurs) composants connectés constitués de tous les points des brins, plus (potentiellement beaucoup) de composants supplémentaires constitués de points de fuite uniques et de ces `` brins de fuite ''. Cependant, ces brins parasites contiennent des points marqués comme «parasites», vous pouvez donc simplement ignorer tout composant contenant un point marqué comme «parasite».

Un danger de cette idée est que vous pouvez avoir une caractéristique où la densité d'un brin diminue progressivement au fur et à mesure que vous vous déplacez le long du brin, jusqu'à ce que la densité soit si faible qu'elle ne soit qu'un ensemble de points parasites. Comme notre critère est «local», il peut ne pas détecter cela et marquer ces points parasites comme faisant partie du brin. Je ne sais pas si ce sera un problème: je suppose que la plupart des points parasites devraient être capturés par le critère, car les changements de densité semblent assez brusques dans votre image.

Si ce problème se produit, vous pouvez essayer une alternative à la simple prise des composants connectés. Pour chaque point X, nous calculons la distance à son plus proche voisin D (X). Nous commençons au point avec un D (X) minimal et effectuons un BFS (ou DFS , l'ordre n'a pas d'importance). Nous ajoutons tout point Y dont D (Y) n'est pas beaucoup plus grand que le D (X) (par un facteur ajustable) avec lequel nous avons commencé. Si nous rencontrons un point Y qui a un D trop grand (Y), nous supprimons le bord (X, Y), marquons Y comme «errant» et agissons comme si nous n'avions jamais visité Y dans notre BFS. S'il est réglé correctement, cela devrait éviter le problème que je décris ci-dessus.

Une autre idée pour résoudre ce problème agit un peu plus localement: vous pourriez faire un BFS et garder une trace du D (X) le plus bas (j'utilise D (X) comme mesure de la densité autour d'un point) rencontré au maximum disons 10 BFS-étapes avant, et si nous rencontrons un Y qui a D (Y) qui est beaucoup plus grand que ce D (X), nous faisons la même chose que l'autre solution (potentielle) que j'ai proposée.

En guise d'avertissement: toutes les idées ci-dessus que j'ai imaginées tout de suite, je ne sais pas vraiment si ce problème particulier a été étudié auparavant, donc je suis peut-être en train de germer un non-sens. Essayez simplement les idées (que ce soit les miennes ou les vôtres) qui vous semblent sensées et découvrez si elles fonctionnent vraiment, puis concentrez-vous uniquement sur leur mise en œuvre efficace.

En utilisant la décomposition modulaire, vous pouvez créer un arbre qui contiendra toutes les particules comme les feuilles et les nœuds supérieurs les regrouperont. Sur la base de cet arbre, vous pouvez définir des mesures qui sont appliquées à chaque nœud de celui-ci, de la racine aux feuilles vers le bas. Vous arrêtez cette traversée vers le bas lorsque les mesures atteignent des seuils définis par l'utilisateur. Une telle mesure peut être la densité de la coque convexe de toutes les particules d'un amas.

Je pense que vous recherchez un algorithme de clustering d'apprentissage automatique.

Cette page de la boîte à outils Python SciKit Learn contient des images qui suggèrent que l' algorithme DBSCAN (Wikipedia) pourrait être ce que vous recherchez. Il semble idéal car son paramètre d'entrée est la taille du quartier, tandis que la plupart des autres algorithmes de clustering veulent le nombre de clusters, que vous ne connaissez pas à l'avance.

"Un algorithme basé sur la densité pour découvrir des clusters dans de grandes bases de données spatiales avec du bruit" Ester, M., HP Kriegel, J. Sander et X. Xu, dans les actes de la 2e conférence internationale sur la découverte des connaissances et l'exploration de données, Portland, OR , AAAI Press, p. 226-231. 1996

J'ai pensé à ce problème. Je ne suis pas un expert en physique, alors soyez indulgent avec moi.

Il semble que ce ne soit pas la distance entre les particules qui compte pour déterminer les mottes. Il s'agit de savoir si les champs de gravité se chevauchent ou non.

Prenez une particule P et déterminez quelles autres particules ont des champs de gravité qui se chevauchent.

Prenez-en un et faites la même chose. Votre objectif n'est pas de trouver toutes les particules dans la touffe mais de trouver ses limites.

Répétez cette opération jusqu'à ce que tous les amas soient trouvés.

Maintenant, revenez en arrière et déterminez la masse des touffes. Vous aurez éliminé les particules parasites, et vous pouvez utiliser les limites des touffes pour trouver la masse.

Je ne sais pas si cela aide, mais c'est tout ce à quoi je pouvais penser.

Vous pouvez, à la fin de chaque pas de temps, convertir les données en un graphique, calculer l'arbre couvrant minimum, puis commencer à supprimer les bords qui dépassent un certain seuil. Cela devrait vous donner des touffes et un moyen facile d'énumérer à travers les particules de chaque touffe.