Astro et programmation

[Leimury]

Salut,

 

A propos de ta base de données, tu l'as segmentée ?

Je pensais à un truc comme une (ou plusieurs) carte par unité de 5°.

[philou31]

Courage bruno...

Juste pour le fun un powerpoint sur le CU7 de Gaia : 1 Terrabytes Petabytes... 10¨12 entités

http://www.pgcon.org/2014/schedule/attachments/334_PGCon2014Gaia.pdf

[jgricourt]

philou31 super intéressant ce case study sur la base de donnée Gaia, voilà un projet sur lequel j'aurai bien aimé bosser ...

['Bruno]

Leimury : j'ai fait des bandes de déclinaison (pas celles du fichier d'origine, d'ailleurs).

 

Philou31 : il va falloir éviter de tout charger en mémoire vive...

[Eric S]

Bonne remarque de Lemury: en théorie on fait d'abord un programme non-optimisé qui marche et qui va servir de référence et ensuite on optimise. La limite pratique, c'est qu'il ne faut pas que la version non-optimisée soit trop lente, genre 10 jours pour afficher la carte .

 

Pour la projection, en photo, on parle habituellement de projection rectilinéaire. C'est simplement la projection (depuis la sphère céleste ici) sur un plan (celui du capteur). La projection, en soit, ne comporte pas de trigo. Et je suis surpris qu'il faille une racine carrée. Historiquement, le calcul de racine carrée a quand même une sale réputation en terme de coût. Je peux imaginer facilement une rotation 3D à partir des (x, y, z) et ensuite sans doute deux divisions (une autre opération à mauvaise réputation mais bien améliorée sur les derniers Intel, ce qui profite aussi à la racine carrée). La rotation 3D, c'est un peu lourd avec 3 x (3 multiplications + 2 additions). Néanmoins, si on reste en champ serré, une des coordonnées va être proche de 1 et on peut l'oublier d'où 2 x (3 multiplications + 2 additions) et plus besoin de division à l'étape suivante.

 

Il faut que j'aille regarder dans mon projet comment j'ai fait ça.

 

philou31: c'est intéressant. Ils utilisent une base de données PostGreSQL mais on ne sait pas s'ils ont des indexations spécifiques pour la répartition sur la voute céleste.

['Bruno]

La racine carrée est nécessaire pour le calcul de la distance entre deux astres.

 

La formule exacte donnant la distance sur la sphère est (de mémoire, donc pas sûr à 100 %) :

cos d = sin δ1 sin δ2 + cos δ1 cos δ2 cos (α2 - α1)

 

Quand la distance est proche de 0, l'arc cosinus est très imprécis, il faut éviter de l'utiliser.

 

Dans ce cas on peut utiliser la formule de Pythagore :

d² = (Δα cos δ)² + (Δδ)²

 

qui correspond d'ailleurs au développement limité à l'ordre 2 de la première formule. La formule de Pythagore est une approximation, mais excellente lorsque la distance est proche de 0. D'où son intérêt pour le dessin à l'oculaire.

 

Le calcul de la racine carrée par l'algorithme de Héron est extrêmement rapide : on double le nombre de décimales exactes à chaque itération. Je ne sais pas si c'est lui qui est utilisé dans l'implémentation, mais si ce n'est pas lui, ç'en est donc un plus rapide (que je ne connais pas). (Tu dis que l'accélération de la division a profité à la racine carrée, donc je suppose que tu as aussi en tête l'algorithme de Héron.)

 

La projection dont tu parles (la sphère céleste sur un capteur) ne serait-elle pas une projection orthogonale ? Pour représenter le ciel elle n'est pas adéquate : trop de déformations.

Modifié 21 septembre 2015 par 'Bruno

[Leimury]

Bonjour,

 

Peut être que tu veux trop mettre les mains dans le cambouis.

Laisses un peu faire l'OS, il est optimisé pour aller tout seul taper sur les copros mathématiques et quelquefois même la carte graphique pour certains calculs.

En tout cas les libs de l'OS sont développées et retouchées sans arrêt depuis 1979 par des générations et des générations d'ingés et de mathématiciens.

 

En fait tu fais quelquefois de la prose sans le savoir

Même si tu bosses en C de base tu n'échapperas pas au C#, pas mal de libs bas niveau sont réécrites pour offrir une adaptation à la machine et quand tu crois tout maitriser en fait tu prends la main sur les optims de l'OS et tu tentes les nationales à la place de l'autoroute.

 

Fais voir un test pour comparer les temps d'un million d'OPs avec ton propre algo de Héron d'un côté et de l'autre en utilisant de bêtes doubles.

 

Evidement n'hésites pas à compiler sans débug et avec toutes les optims de vitesse (aucune optim de taille).

Vérifies également les optims maths

Et bien sûr avec un petit tirage au sort des valeurs pour éviter de ressasser 1000000000000 fois le même résultat stocké.

 

Struct Member align : 4 ou 8

SecurityCheck et Basic runtime check tu dégages.

Floating point model : Tu as le choix entre précise, fast ou strict.

 

En débug, sans aucune optim, il est possible qu'un algo maison dépote par rapport au truc tout bête et facile à lire.

 

Compilé en final ou release avec toutes les optims qui vont bien c'est fou comme on se rend compte que la facilité de lecture peut aussi avoir du bon.

 

En fait plus tu restes dans des schémas habituels et plus l'optimisation tend à marcher.

 

Je parie que tu compiles en ligne de commande...

Vu que tu développes en C sous Windows il pourrait être utile de vérifier toutes les options de compil que tu trouveras facilement dans les panneaux d'options de VisualStudio Express : Ca peut donner des idées.

 

Les mots clefs du type Register sont à éviter.

Par contre les optims à base de zones ou variables statiques sont toujours aussi efficaces quand c'est utilisable.

Evidement moins ta fonction crée de variables et moins tu as d'overhead.

 

Fais également gaffe au multithreading.

Si tu programmes en monothread et que tu compiles en multi-threading, tu utiliseras des libs plus lourdes et la compil produit quelque chose de plus lent.

Si tu as du /MT ou /MTd (MultiThread et MultiThread Debug) dans tes options c'est pas bon si tu bosses en mono thread.

 

Bon ciel

Modifié 21 septembre 2015 par Leimury

[Fred_76]

Bonjour

 

As tu regardé les sources de Stellarium ?

 

Ils gèrent eux aussi des gros catalogues d'étoiles et l'affichage est fluide.

 

À+

 

Fred

['Bruno]

Leimury : tu n'as pas compris mon message. Je n'ai jamais dit que j'allais écrire moi même l'algorithme de Héron, j'ai dit que le calcul de la racine carrée est rapide puisqu'il existe un algorithme rapide pour ça ─ soit ça se fait avec Héron, soit se fait autrement et alors c'est encore plus rapide.

 

Bref, pour récapituler : j'ai écrit un programme qui génère un affichage fluide (Fred_76 : je n'ai pas pensé à regarder Stellarium parce que je n'en ai pas ressenti le besoin et que je n'y ai pas pensé) qui n'est pas le programme définitif, mais le programme m'ayant servi à tester diverses méthodes, notamment pour optimiser (Leimury : inutile de me conseiller d'écrire le programme puis de regarder ce qui prend du temps pour optimiser, c'est précisément ce que j'ai fait). C'est une fois ce premier programme terminé que j'ai eu l'idée d'en parler ici pour donner suite à mon vieux message. J'ai encore deux ou trois optimisations à faire mais rien n'est sûr (refaire des bandes de déclinaisons plus petites ; stocker les cosinus et les sinus des coordonnées ? ça va augmenter la taille des tableaux...) [Et je n'ai pas encore fait l'affichage des nébuleuses parce que je ne sais pas comment m'y prendre pour l'instant.]

 

---

Finalement, par curiosité, je viens d'aller voir les sources de Stellarium. Plus de 100 Mo ! Des centaines de fichiers .cpp (en plus je ne connais pas le C++...) et une énorme hiérarchie de répertoires. Bref, c'est du lourd, ça me prendrait un temps fou pour comprendre comment c'est organisé...

Modifié 21 septembre 2015 par 'Bruno

[Eric S]

Je viens de regarder dans mon projet. Le code n'est pas optimisé mais néanmoins ni racine carrée, ni arcsin/cos.

 

Je le mets ici si d'autres sont intéressés:

 

 Scale := sqrt(sqr(Bitmap.Width)+sqr(Bitmap.Height)) / (2.0 * tan(FOV2));

 RotX := TRotationMatrix.Create;
 RotY := TRotationMatrix.Create;
 RotZ := TRotationMatrix.Create;
 try
   // euler transformation, 
   RotZ.Rz(-RA);
   RotY.Ry(Decl);
   RotX.Rx(North);

   for I:=low(Stars) to high(Stars) do
   begin
     if Stars[i].MagV<MagSat then
       Gray := 255
     else
       Gray := round(255.0*power(100,-Stars[i].MagV/5.0)/power(100,-MagSat/5.0));
     Bitmap.Canvas.Brush.Color:=RGBToColor(Gray, Gray, Gray);
     Bitmap.Canvas.Pen.Color:=Bitmap.Canvas.Brush.Color;

     SinCos(Stars[i].Decl, SinDecl, CosDecl);
     SinCos(Stars[i].RA, SinRA, CosRA);

     Vect.coef[0] :=  cosRA * cosDecl;
     Vect.coef[1] := -sinRA * cosDecl;
     Vect.coef[2] := sinDecl;


     Vect := RotZ.Rotate(Vect);
     Vect := RotY.Rotate(Vect);
     Vect := RotX.Rotate(Vect);

     // projection on the screen
     Inv := Scale/Vect.coef[0];
     Vect.coef[1] := Vect.coef[1] * Inv;
     Vect.coef[2] := Vect.coef[2] * Inv;

     Bitmap.Canvas.EllipseC(round(Vect.coef[1] + 0.5 * Bitmap.Width),
                            round(0.5 * Bitmap.Height - Vect.coef[2]),
                            round(Blur_Radius), round(Blur_Radius));
   end;

 

Je fais trois rotations 3D successives qu'on pourrait optimiser en amont en une seule.

 

Je recalcule les sin/cos dec/ad puis les x, y et z à chaque fois. Je pourrais stocker les x, y et z directement.

 

Pour la projection, c'est la projection rectilinéaire, cad depuis le centre de la sphère. Elle a tendance à agrandir loin du point de contact entre la sphère et le plan.

 

La projection orthogonale projette perpendiculairement au plan et compresse les éléments lointains.

 

Enfin les cartographes utilisent la projection stéréographique qui projette sur le plan depuis le point de la sphère à l'opposé du contact entre le plan et la sphère. Elle est intermédiaire entre les deux premières et agrandi légèrement les bords (elle conserve aussi les angles).

 

Pour un angle de vue de 20°, la différence entre les deux projections extrêmes est de 1,5% seulement. Donc on peut facilement faire des approximations avec des champs assez larges.

 

Dans mon code, ça donnerait Vect.Coef[0]~1 donc pas besoin de le calculer et

 

     
// projection on the screen
     Vect.coef[1] := Vect.coef[1] * Scale;
     Vect.coef[2] := Vect.coef[2] * Scale;

[Leimury]

Leimury : tu n'as pas compris mon message. Je n'ai jamais dit que j'allais écrire moi même l'algorithme de Héron' date=' j'ai dit que le calcul de la racine carrée est rapide puisqu'il existe un algorithme rapide pour ça ─ soit ça se fait avec Héron, soit se fait autrement et alors c'est encore plus rapide.[/quote']

 

Désolé, j'avais pas compris.

 

refaire des bandes de déclinaisons plus petites ; stocker les cosinus et les sinus des coordonnées ? ça va augmenter la taille des tableaux...)

 

N'aies pas trop de scrupules avec la tailles des données.

N'oublies pas que 2Go c'est qu'une toute petite clef de rien du tout.

 

En tout cas bonne chance pour la suite

[Fred_76]

Pour trier tes étoiles à afficher, je pensais à une hierarchie par magnitude et proximité.

 

Chaque étoile de forte magnitude (il n'y en a pas beaucoup) possède les références des étoiles de moindre magnitude qui l'entourent. Puis chacune de moindre magnitude celles d'encore moindre magnitude. Cela sur 3, 4 ou 5 niveaux voire plus.

 

La sélection commence donc par une recherche sur le premier niveau, ce qui va assez vite car il y a peu de membres à tester. Ensuite pour le second niveau, et ainsi de suite. A chaque fois le nombre de membres à chercher est relativement faible d'autant plus qu'il ne faut se concentrer que sur les membres proches de la bordure du champ à afficher. Ceux qui sont complètement dans le champ seront intégralement affichés sans autre calcul.

 

Je ne sais pas si tu vois à quoi je pense !

 

La préparation des données est un programme à part entière...

['Bruno]

Éric S : très intéressant, je vais examiner ça à tête reposée !

 

A priori je choisissais la projection pour diminuer les déformations, mais si celles-ci sont peu importantes, effectivement il faut que je pense à choisir les projections qui diminuent les calculs. (Remarque : la "projection" qui utilise des racines carrées n'utilise aucune fonction trigo, pas même directe. Je pense qu'elle reste la plus rapide.)

 

Fred_76 : excellente idée ! Quand bien même préparer les fichiers qui utilisent cette hiérarchie serait long et compliqué, ce serait fait une fois pour toutes. (Pour l'instant je stocke de la plus brillante à la plus faible uniquement pour que les étoiles faibles s'affichent par dessus les étoiles brillantes au lieu d'être effacées par elles, mais ton idée est compatible.) C'est une méthode qui me paraît compliquée, mais je suis persuadé de l'intérêt de préparer les données au mieux. Maintenant que je suis dans la phase où je veux tout refaire proprement, ce genre de réflexion m'intéresse !

Modifié 21 septembre 2015 par 'Bruno

[Fred_76]

Ca m'a servi ces artifices de programmation, il y a une 20aine d'années, quand on m'a demandé d'optimiser un logiciel sous IBM 36, en Fortran, qui prenait 3 à 4 heures de calcul pour faire 10 simulations. Après optimisation, il ne prenait plus que 1 à 2 minutes pour en faire 1 000 ! Et quand on a compilé ça sur un PC, il faisait 10 000 en moins de 30 s.

 

L'artifice des précalculs des sin x, cos x, sin 2x, cos 2x... a été génial, mais c'est surtout le contrôle des précisions et des boucles qui a été crucial, pour ne pas perdre du temps à calculer des choses inutilement précises, ou tout simplement inutiles. Il faut donc se donner un objectif et s'y tenir sans vouloir faire de la surqualité.

 

Par contre le nombre de lignes de code du programme a été démultiplié...

[Eric S]

Je me permets d'insister parce que mon logiciel n'est pas optimisé et rame quand je le fais tourner sur un PDA (sorti en 2008 certes...) . Même si j'ai abandonné le développement de mon logiciel, ça pourrait toujours servir plus tard pour une autre application ou quand on va se prendre le gros catalogue produit par Gaïa .

 

Le calcul le plus simple (sans la rotation de champ), c'est:

X = (AD - AD_0) / cos (Dec_0)

Y = Dec - Dec_0

Pas de trigo, pas racine carrée...

... mais on sent bien qu'il y a un gros problème vers le pôle avec cos(Dec_0) qui tend vers 0 par exemple.

 

Néanmoins:

 

Mais aussi parce que je calcule la distance angulaire avec :

cos(distance) = sin b x sin b0 + ( cos b x cos b0 x cos(l - l0))

[...]

 

J'imagine ce que tu peux faire avec la distance angulaire' date=' c'est-à-dire dessiner ton étoile depuis le centre de la carte. Mais alors comment calcules-tu la direction de l'étoile?

 

 

Sinon, pour les optimisations, tu peux juste stocker dans ton fichier les AD/Dec et calculer les sin/cos et/ou x/y/z la première fois que tu les utilises. Ça limite la taille sur le disque dur et le temps de lecture. Le premier affichage va être un peu lent mais moins que de tout calculer à l'ouverture programme. Par contre, l'empreinte en mémoire vive va rester la même (légèrement augmenter vu qu'il faut un booléen pour savoir si le calcul a déjà été fait).

 

Quant à charger tout le fichier en mémoire, ça se discute aussi. Les OS actuels vont de toute façon mettre les morceaux du fichier déjà lus en mémoire cache. C'est seulement quand tu changes de zone qu'il risque d'y avoir des ralentissements. Par exemple, avec [i']astrometry.net[/i] en local, la première résolution est un peu lente mais ensuite les fichiers d'index étant déjà en mémoire, ça va plus vite sans que rien dans le code du logiciel ne force à garder les données en mémoire.

['Bruno]

Le calcul le plus simple (sans la rotation de champ), c'est:

X = (AD - AD_0) / cos (Dec_0)

Y = Dec - Dec_0

Pas de trigo, pas racine carrée...

C'est celui que j'utilise pour la vue à l'oculaire, c'est une projection rectangulaire, et en effet il n'y a pas de racines carrées (je ne sais plus pourquoi j'ai parlé de racines carrées, tiens, ah oui : comme approximation du calcul de la distance angulaire, mais ici ça ne sert pas). Mais attention : ça affiche une carte avec le nord en haut. Moi je veux le zénith en haut. Je dois donc ajouter une rotation.

 

... mais on sent bien qu'il y a un gros problème vers le pôle avec cos(Dec_0) qui tend vers 0 par exemple.

Près du pôle il faut une projection polaire. Même si le pôle n'est pas dans le champ mais un peu à côté, les déformations de la projection rectangulaire deviennent trop importantes. Donc il faut calculer une distance (avec Pythagore, pas besoin de la formule exacte) et un angle.

 

J'imagine ce que tu peux faire avec la distance angulaire, c'est-à-dire dessiner ton étoile depuis le centre de la carte.

C'est comme ça que je procède pour la vue à l'oeil nu et la vue au chercheur.

 

Mais alors comment calcules-tu la direction de l'étoile?

En calculant l'angle de position ! (C'est analogue aux coordonnées polaires : on a une distance et un angle à calculer.)

 

Sinon, pour les optimisations, tu peux juste stocker dans ton fichier les AD/Dec et calculer les sin/cos et/ou x/y/z la première fois que tu les utilises. Ça limite la taille sur le disque dur et le temps de lecture.

La taille qui me préoccupe, c'est plus celle de la mémoire vive. Cela dit je pense que tu as raison, du moins si calculer les sinus et cosinus n'est pas beaucoup plus lent que lire sur le disque dur, ce que je soupçonne (à essayer).

 

Quant à charger tout le fichier en mémoire, ça se discute aussi.

C'est quand même beaucoup plus simple à gérer...